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ABSTRACT 

This  report  lays  out  the  mathematical  framework  and  reasoning  involved  in 
addressing  the  question  of  how  to  produce  sophisticated  false  targets  in  both 
range  and  Doppler  to  jam  a  radar.  The  explanations  here  are  given  from  first 
principles  in  order  to  clarify  the  approach  used  in  writing  software  to  address 
the  task.  They  constitute  standard  radar  theory  seen  through  a  physicist’s 
eyes,  and  are  recorded  as  background  to  the  approach  taken  in  addressing  the 
jamming  task  in  a  later  report. 

We  discuss  in  detail  how  a  radar  generates  a  range-Doppler  plot,  using 
a  set  of  parameters  that  describe  the  outgoing  radar  signal  and  the  relevant 
characteristics  of  the  target  whose  impulse  response  is  known.  We  also  ex¬ 
plain  necessary  concepts  from  a  mathematical  physicist’s  viewpoint:  bounds 
on  pulse  parameters,  correlation/convolution,  the  theory  of  Hilbert  transforms, 
the  relevant  Fourier  analysis,  and  general  concepts  of  radar  signal  processing 
such  as  ambiguity  functions  and  the  maximum  detectable  range  of  a  target. 
This  entire  discussion  is  aimed  at  indicating  the  approach  and  philosophy  used 
to  solve  the  various  problems  encountered  while  working  on  the  task. 
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How  to  Create  and  Manipulate  Radar  Range— Doppler  Plots 

Executive  Summary 

This  report  lays  out  the  mathematical  framework  and  reasoning  involved  in  addressing 
the  question  of  how  to  produce  sophisticated  false  targets  in  both  range  and  Doppler  to 
jam  a  radar.  The  explanations  here  are  given  from  first  principles  in  order  to  clarify  the 
approach  used  in  writing  software  to  address  the  task.  They  constitute  standard  radar 
theory  seen  through  a  physicist’s  eyes,  and  are  recorded  as  background  to  the  approach 
taken  in  addressing  the  jamming  task  in  a  later  report,  currently  in  publication. 

We  first  derive  the  relevant  radar  equations  that  allow  the  radar  scenario  to  be  modelled 
numerically.  An  in-depth  discussion  of  correlation  follows,  which  is  central  to  the  operation 
of  a  matched  filter.  Next  we  explain  how  to  model  the  returned  signal  from  knowledge  of 
the  target’s  impulse  response.  Following  a  discussion  of  Doppler  windowing  is  an  analysis 
of  what  constitutes  valid  and  useful  pulse  and  processing  parameters,  such  as  sampling 
interval  and  bandwidth,  with  an  explanation  of  blind  speed  and  velocity  aliasing.  We 
then  explain  the  mathematics  of  how  a  jammer  can  add  time-dependent  phases  to  a  signal 
with  the  aim  of  manipulating  that  signal’s  Doppler  content,  together  with  the  simpler 
manipulation  of  its  range  content  carried  out  by  adding  delay  to  the  signal. 

Appendices  discuss  the  Hilbert  transform  and  the  relation  of  correlation  to  convolution, 
including  first-principles  examples  of  these.  Calculation  of  signal-to-noise  ratios  and  how 
to  combine  the  cross  sections  of  many  scatterers  is  also  included.  The  report  closes  with 
an  explanation  of  the  use  of  the  ambiguity  function.  Explanations  of  textbook  concepts  of 
radar  signal  processing  have  been  given  in  this  report  as  a  way  of  indicating  the  approach 
and  philosophy  used  to  solve  the  various  problems  encountered  while  working  on  the 
jamming  task. 
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1  Introduction 

This  report  lays  out  the  mathematics  involved  in  simulating  the  jamming  of  a  radar  to 
induce  it  to  generate  false  targets  having  structure  in  both  range  and  Doppler.  We  begin 
with  the  basic  concepts  of  radar  signal  processing,  then  build  on  them,  step  by  step,  to 
arrive  at  the  relevant  false-target  mathematics.  This  mathematics  can  then  be  used  to 
simulate  numerically  the  interaction  of  a  radar  with  a  target  composed  of  perhaps  many 
thousands  of  scatterers  of  known  cross  section.  By  simulating  how  the  radar  interacts,  we 
can  then  also  simulate  jamming  it  to  produce  false  targets. 

We  start  by  setting  up  the  relevant  radar  equations  that  allow  the  radar  scenario  to 
be  modelled  numerically.  An  in-depth  discussion  of  correlation  follows,  which  is  central 
to  the  operation  of  a  matched  filter.  After  a  discussion  of  Doppler  windowing  and  why  it 
is  generally  used,  is  an  analysis  of  what  constitutes  valid  and  useful  pulse  and  processing 
parameters  such  as  sampling  interval  and  bandwidth,  with  an  explanation  of  blind  speed 
and  velocity  aliasing.  We  then  explain  how  time-dependent  phases  can  be  added  to  a 
“returned”  signal  sent  by  a  jammer  to  a  radar  with  the  aim  of  manipulating  that  signal’s 
target- velocity  content,  together  with  the  simpler  manipulation  of  its  range  content  which 
can  be  done  by  adding  delay  to  the  signal.  Finally,  appendices  discuss  the  Hilbert  trans¬ 
form,  the  relation  of  correlation  to  convolution  (with  examples  of  how  correlation  and 
convolution  are  implemented  from  first  principles),  signal- to- noise  theory,  and  the  use  of 
the  ambiguity  function. 


2  The  Mathematics  of  Emitted 
and  Received  Signals 

The  pages  that  follow  show  each  step  of  assembling  the  mathematics  that  describes  a  radar 
interaction.  We  begin  with  a  discussion  of  the  use  of  complex  numbers  to  represent  a  real 
signal. 

The  electromagnetic  field  that  comprises  a  propagating  radar  signal  is  very  complicated 
both  in  space  and  in  time.  But  there  is  no  requirement  to  analyse  it  exactly:  we  need  only 
consider  its  far-held  component,  and  of  this,  it’s  sufficient  to  work  with  the  electric  held 
only,  since  that  is  always  related  to  the  magnetic  held  in  a  well-defined  way.  At  any  given 
point,  the  electric  held  of  a  radar  signal  can  be  represented  by  a  unit-length  vector  (the 
held’s  polarisation)  which  is  multiplied  by  a  number  that  varies  sinusoidally  in  time  with 
frequency  /.  We  will  assume  the  vector  is  given  (i.e.  the  polarisation  is  known),  and  will 
represent  the  held  by  the  sinusoidally  varying  number  alone. 

Sinusoids  are  particularly  useful  for  signal  processing,  being  a  special  case  of  signals  with  an 
exponential  form.  An  exponential  signal  is  unique  in  that  its  form  is  preserved  when  it  passes 
through  a  linear  time-shift-invariant  system.  A  linear  system  can  be  treated  as  processing  a  sum  of 
signals  component-wise  and  adding  the  results;  a  time-shift-invariant  system  is  one  that  processes 
identically  today  as  it  did  yesterday.  Most  of  the  important  systems  in  the  signal  processing  world 
are  linear  time-shift  invariant.  The  usefulness  of  the  exponential  signal  here  embodies  the  fact 
that  exponential  signals  are  the  eigenfunctions  of  linear  time-shift-invariant  systems. 

The  information  content  of  a  signal  is  almost  always  placed  as  a  modulation  onto  a  high- 
frequency  carrier  wave.  The  use  of  this  carrier  is  only  partly  dictated  by  available  wave¬ 
generating  technology  and  transparency  of  Earth’s  atmosphere  to  different  wavelengths. 
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One  main  reason  for  transmitting  information  on  a  high-frequency  carrier  follows  from 
the  fact  that  waves  of  length  A  emitted  coherently1  from  an  aperture  of  diameter  a  will 
diffract  to  form  a  beam  of  angular  width  approximately  A/a.  Only  by  using  relatively 
high  frequencies  (meaning  small  A)  will  a  beam’s  width  be  small  enough  that  it  will 
deposit  a  good  fraction  of  its  energy  on  a  distant  target.  Another  reason  for  using  high 
frequencies  follows  from  the  fact  that  the  radio  waves  travel  at  a  set  speed.  To  transmit 
more  information,  we  can’t  hope  to  send  the  same  wave  at  a  higher  speed.  Instead,  we 
must  place  a  higher  density  of  information  into  the  wave,  by  giving  the  modulation  more 
structure:  and  this  we  do  by  shortening  the  duration  of  the  wave’s  modulation  changes  that 
encode  the  information.  Thus  the  wave’s  shape  must  be  made  to  change  more  rapidly, 
and  Fourier  analysis  tells  us  that  such  rapid  changes  require  the  signal  to  be  a  sum  of 
monochromatic  waves  over  a  large  spread  of  frequencies.  (That  is,  high  data  rates  are 
built  from  large  bandwidth .)  The  bandwidth  is  spread  around  the  carrier  frequency,  and 
if  the  frequencies  making  up  a  large  bandwidth  are  all  to  be  high  enough  to  ensure  that 
the  beam  has  the  required  low  angular  width,  then  the  carrier  frequency  must  be  high. 

Yet  another  reason  for  using  a  high  carrier  frequency  involves  the  range-velocity  trade¬ 
off  examined  on  page  34:  higher  carrier  frequencies  enable  the  radar  to  measure  both  larger 
ranges  and  higher  radial  velocities  in  a  single  measurement.  But  these  higher  frequencies 
also  produce  more  “blind  speeds”,  discussed  on  page  35. 

Suppose  then  that  we  have  formed  a  signal  by  modulating  both  the  amplitude  and 
phase  of  a  carrier  wave  of  angular  frequency  cuq  =  27r/o.  We  have  given  the  signal  a  large 
bandwidth,  but  in  practice  this  bandwidth  is  still  much  less  than  the  carrier  frequency.  In 
that  case,  the  amplitude  and  phase  do  not  vary  so  drastically  as  to  alter  the  signal  radically 
from  a  recognisable  sinusoid;  in  other  words,  the  signal  is  confined  to  a  relatively  narrow 
frequency  interval  around  (and  compared  to)  its  carrier.  This  type  of  signal  is  often  referred 
to  as  narrow  band ;  even  though  it  might  have  a  large  bandwidth,  it  is  “relatively”  narrow 
band  compared  to  the  carrier.  (In  contrast,  a  signal  that  is  not  narrow  band  is  not  usefully 
described  by  the  “analytic  signal”  analysis  that  follows.)  We  can  represent  a  narrow-band 
signal  by  a  sine  or  cosine;  the  choice  is  arbitrary  since  the  two  functions  differ  only  by  a 
constant  phase.  Suppose  we  choose  a  sine,  writing  the  signal  as  A(t)  sin[ca0i  +  <(>(i)]  where 
A{t)  and  if>(t)  are  the  signal’s  instantaneous  amplitude  and  phase  respectively.  (We’ll  drop 
the  time  dependence  of  A  and  <)>  for  brevity  in  the  equations  that  follow.)  The  task  of 
the  radar  receiver  is  to  determine  A  and  <j)  from  one  measurement  to  the  next.  Note  that 
“phase”  is  sometimes  taken  to  mean  </>(f )  and  sometimes  u>o t  +  4>{t),  but  there  should  be 
no  ambiguity  in  the  discussion  that  follows. 

Focus  on  determining  i f>,  since  we’ll  see  that  the  amplitude  A  emerges  along  with  cj). 
The  phase  is  an  angle,  and  an  angle  always  needs  two  pieces  of  information  to  be  fully 
determined.  Certainly  any  angle  is  just  one  number,  but  it’s  the  ratio  of  two  quantities:  the 
arc  length  that  the  angle  scribes  on  a  given  circle,  and  the  radius  of  that  circle.  Ordinarily 
the  radius  is  set  to  have  length  one,  in  which  case  the  angle  is  a  single  number,  the  length 
of  the  arc — but  this  is  a  convention.  Another  way  of  realising  that  an  angle  needs  two 

xBy  coherent  waves  is  meant  a  set  of  wave  fronts  with  constant  wavelength,  direction,  and  relative 
phase.  Electromagnetic  waves  produced  by  radar  transmitters  and  lasers  are  coherent;  electromagnetic 
waves  produced  by  a  household  torch  are  not.  The  waves  flowing  from  a  torch  are  coherent  over  tiny 
time  scales  and  tiny  distances,  but  this  coherence  is  constantly  changing  and  so  is  somewhat  trivial.  Real 
coherence  is  all  about  the  large-scale  effects  of  wave  interference,  so  it  requires  “tidy”  waves  at  large  scales 
and  long  periods  of  time. 
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pieces  of  information  to  be  fully  determined  is  to  represent  the  angle  by  a  vector  of  any 
length  r  in  the  xy  plane;  the  vector’s  two  components,  r  cos  (j)  and  r  sin  cj),  are  the  required 
two  pieces  of  information  that  fully  determine  (j>.  In  fact  the  length  r  is  redundant,  since 
it  can  be  extracted  from  r  cos  <f>  and  r  sin  (ft  by  summing  the  squares  of  these  two  numbers. 
So  we  can  ignore  it,  and  our  two  pieces  of  information  that  will  fully  determine  the  phase 
of  the  signal  will  be  cos  (j)  and  sin  cf>.  In  principle,  these  can  be  extracted  in  the  following 
way  when  representing  the  signal  by  a  sine  or  by  a  cosine.  (More  discussion  of  this  in  a 
radar  hardware  context  can  be  found  in  [1].) 


Representing  the  Signal  by  a  Sine:  Write  the  signal  as  Asin(cJot  +  4>)  and  ex¬ 
tract  cos  (j)  and  sin  f>  by  mixing  (multiplying)  the  signal  with  reference  signals  2sinix>ot 
and  2  cos^ot.  In  equations  (2.2)-(2.8)  that  follow,  we  will  employ  the  following  standard 
trigonometric  identities: 

—2  sin  a  sin  (3  =  cos(a  +  (5)  —  cos(a  —  (3) ,  2  sin  a  cos  (3  =  sin(ci!  +  /3)  +  sin(ci!  —  /3 ) , 

2  cos  a  cos  (3  =  cos(a  +  f3)  +  cos(a  —  /3) ,  cos(a  +  /3)  =  cos  a  cos  (3  —  since  sin/3  .  (2.1) 

First,  mix  the  signal  with  2sinwot: 

A  sin(w0t  +  <f>)  2  sin  u0t  =  A  cos  <f>  —  A  cos(2w0f  +  <fi)  .  (2-2) 

i _ _j  i _ i 

“I  data”  filtered  out 

The  second  term  on  the  right-hand  side  of  (2.2)  has  a  high  frequency  and  so  is  removed 
by  a  low-pass  filter.  The  first  term  is  our  first  piece  of  information  that  determines  </>. 
Because  it  was  produced  by  mixing  our  representative  signal  (a  sine)  with  an  oscillation 
described  in  the  same  way  (i.e.  also  a  sine),  it’s  called  the  in-phase  data  (“I  data”)  of  the 
signal. 

Now  mix  the  signal  with  2  cos  uj0t: 

A  sin(woi  +  (j>)  2  cos  u0t  =  A  sin  <j>  +  A  sin(2o;ot  +  <j>)  ■  (2-3) 

i - j  i _ i 

“Q  data”  filtered  out 

Again  the  second  term  is  filtered  out,  and  the  first  term  is  called  the  signal’s  quadra- 
phase  data 2  (“Q  data”).  Traditionally,  the  I  and  Q  data  are  combined  into  a  phasor:  a 
vector  with  xy  coordinates  ( I,Q )  =  A(cos  4>,  sin  cff).  The  amplitude  A(t)  of  the  signal  is 
the  time-varying  length  of  this  phasor. 

The  above  mixing  process  is  reversed  to  convert  the  I/Q  data  back  to  a  modulated 
signal  ready  to  be  sent  by  a  radar  transmitter.  The  I  data  in  (2.2)  resulted  from  mixing 
the  signal  with  a  sine,  so  now  mix  the  I  with  a  sine;  the  Q  data  in  (2.3)  resulted  from 
mixing  the  signal  with  a  cosine,  so  now  mix  the  Q  with  a  cosine.  Then  add  the  results: 

A  cos  sin  co0t  +  A  sin  4>  cos  uj0t  =  A  sin(u;of  +  4>) .  (2-4) 

i _ i  i _ i 

I  data  Q  data 

The  resulting  signal  Asin(u;ot  +  <j>)  can  then  be  sent  out  by  the  transmitter. 

2The  use  of  “quad”  here,  as  in  quadrant  or  square,  relates  to  the  idea  of  a  square’s  right-angled  corners: 
these  mimic  the  right  angle  between  two  phasors  that  represent  a  sine  and  a  cosine. 
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Representing  the  Signal  by  a  Cosine:  If  we  choose  instead  to  write  the  signal 
as  Acos(ujQt  +  p)  (where  this  p  differs  from  that  used  in  the  sine  representation  above), 
very  little  changes  in  the  above  analysis.  Again  we  mix  the  signal  with  2  sin  c o^t  and  2  cos  u^t. 
This  time  begin  with  2  cos  u0t: 

A  cos (uj0t  +  p)  2  cos  u0t  =  A  cos  p  +  A  cos(2u;o t  +  p)  .  (2-5) 

i - 1  i _ _ _ i 

I  data  filtered  out 

The  second  term  is  filtered  out,  and  because  the  first  term  was  produced  by  mixing  our 
representative  signal  (a  cosine)  with  an  oscillation  described  in  the  same  way  (i.e.,  a  cosine), 
it’s  called  the  I  data  of  the  signal,  just  as  for  the  sine  case  above.  Now  mix  the  signal 
with  2  sin  Wot: 

A  cos(w0t  +  (p)  2  sin  co0t  =  —  A  sin  p  +  A  sin(2w0£  +  p)  .  (2-6) 

i - 1  i _ i 

Q  data  filtered  out 

Filter  the  second  term  out  and,  as  before,  call  the  first  term  the  signal’s  Q  data.  Again 
the  I  and  Q  data  can  be  combined  into  a  phasor,  but  this  time  a  minus  sign  is  included: 
(I,—Q)  =  A(cos  p,  sin  p).  We  see  that  the  resulting  phasor  is  identical  to  that  formed 
when  representing  the  signal  by  a  sine  (that  is,  if  we  include  the  minus  sign). 

Again  the  I/Q  data  can  be  converted  to  a  signal  ready  for  transmitting.  The  I  data 
in  (2.5)  resulted  from  mixing  the  signal  with  a  cosine,  so  now  mix  the  I  with  a  cosine;  the 
Q  data  in  (2.6)  resulted  from  mixing  the  signal  with  a  sine,  so  now  mix  the  Q  with  a  sine. 
Then  add  the  results: 

A  cos  p  cos  u0t  —  Asin^  sinwot  =  Acos(u;o£  +  0) .  (2-7) 

i _ i  i _ i 

I  data  Q  data 

The  resulting  signal  Acos(wgt  +  p)  is  now  sent  out  by  the  transmitter. 

The  above  mixing  procedure  is  all  about  our  choosing  to  represent  a  sinusoid  as  either 
a  sine  or  a  cosine.  We  must  choose  one  or  the  other,  and  the  choice  corresponds  to  our 
deciding  what  to  label  as  the  I  data  and  what  to  label  as  the  Q  data  after  the  mixing 
process.  “I  data”  refers  to  the  output  of  mixing  the  signal  with  the  sinusoid  that  will 
represent  the  signal;  i.e.,  to  the  output  of  mixing  the  signal  with  2sinwot  if  we  choose 
to  represent  the  signal  by  a  sine,  and  2cosco0t  if  we  choose  to  represent  the  signal  by  a 
cosine.  The  I  and  Q  numbers  then  describe  the  signal  completely. 

Refinements  from  Superhet  Receivers  In  practice,  most  radar  receivers  are  of  the 
superheterodyne  type;  instead  of  removing  the  carrier  in  one  step,  superheterodynes  remove 
it  in  a  series  of  smaller  steps  as  a  way  of  improving  the  final  signal-to-noise  ratio.  But  the 
procedure  at  each  step  is  identical  to  that  described  above.  For  example,  begin  with  a 
signal  Asin(w0t  +  p)  and  mix  it  with  a  reference  wave  that  incorporates  an  intermediate 
frequency  wjF: 

Asin(w0t  +  i p)  2cos[(w0  —  uqF)£]  =  Asin(wIFt  +  p)  +  Asin[(2w0  —  u;IF)t  +  p]  .  (2.8) 

The  second  term  on  the  right-hand  side  of  (2.8)  is  now  filtered  out,  leaving  the  first 
term  which  is  just  like  the  original  signal  but  with  a  new,  lower- frequency  carrier  uqp. 
The  mixing/filtering  process  is  now  applied  to  this  term  (using  a  still  lower  intermediate 
frequency)  in  one  or  more  steps  until  we  have  removed  the  carrier  entirely.  Along  the  way 
several  stages  of  amplification  will  also  have  been  applied. 
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2.1  Using  Complex  Numbers  to  Represent  I/Q  Data 


The  numbers  A  and  fi — really  A(t),fi(t ) — give  a  phasor  A(cosfi,  sinfi)  whose  relatively 
small  changes  over  time  embody  the  modulation  of  the  carrier  wave.  A  phasor  of  the  actual 
radio-frequency  signal  can  be  constructed  by  making  the  A-(f>  phasor  spin  at  constant 
angular  speed  uj0.  To  see  why,  simply  rotate  the  A-fi  phasor  by  angle  co0t  using  a  matrix 
multiplication: 


spinning  phasor 


cos  u0t 

—  sin  uj0t 

A  cos  fi 

A  cos  (uj0t  +  fi) 

sin  u0t 

cos  u0t 

A  sinfi 

A  sin(a;o^  +  fi) 

_ i  i _ i  i _ i 

rotate  vector  by  cj0£  vector  to  be  rotated  vector 

rotated 


(2.9) 


so  that  the  radio-frequency  signal  is  either  the  x-  or  the  y-component  of  this  rapidly 
spinning  phasor. 

Any  phasor  is  traditionally  a  vector  [  £  ] ,  but  converting  it  to  the  complex  number  a  +  ib 
makes  it  arguably  easier  to  manipulate.  For  example,  in  the  vector  picture,  to  increase 
the  phase  by  angle  #  we  rotate  the  phasor  vector  using  a  matrix  multiplication: 


a 

cos  9 

—  sin# 

a 

b 

rotated  by  6 

sin# 

cos  # 

b 

(2.10) 


in  contrast,  the  complex-number  version  of  (2.10)  employs  “everyday”  multiplication: 


(a  +  ib) rotated  by  6  =  (cos  9  +  i  sin  9)  (a  +  ib) .  (2.11) 

This  complex  notation  is  compact  but  not  mystical.  Equation  (2.11)  is  a  “rendering  down” 
of  (2.10)  in  which  the  seemingly  redundant  repetition  of  the  sine  and  cosine  in  the  rotation 
matrix  has  been  removed,  at  the  small  cost  of  introducing  an  “i”.3 4 * *  The  term  cos  9  +  isin  9 
is  written  as  eid  by  defining  the  exponential  of  a  complex  number  to  have  the  same  Taylor- 
series  form  as  the  exponential  of  a  real  number.  Complex  notation  converts  (2.9)  to 

“circling”  complex  number  =  j  _  ,ie1^ ,  =  i  Ael^°t+^  i ,  (2.12) 

rotate  complex  complex  number  rotated  complex 
number  by  ujQt  to  be  rotated  number 

which  can  (arguably)  be  thought  of  as  simplifying  the  mathematics  of  rotating  phasors: 
compare  the  elaborate  layout  of  (2.9)  with  the  simplicity  of  (2.12). 

The  process  of  extracting  the  I/Q  content  at  time  t  from  a  rotating  phasor  can  now  be 
viewed  as  multiplying  it  by  which  removes  its  rapidly  rotating  carrier  content  (or 

equivalently,  switches  to  a  rotating  frame  in  phasor-space!).  Likewise,  converting  I/Q  data 
to  a  signal  ready  to  send  out  at  time  t  can  be  viewed  as  multiplying  its  phasor  by  e1 *^0*. 

3Here  we  have  rendered  down  a  two-dimensional  rotation  matrix.  When  the  same  idea  is  applied  to 
three-dimensional  rotations,  what  results  is  quaternions,  whose  i,j,k  seem  to  attract  mystical  descriptions 
by  some  in  the  computer  graphics  community.  But  again,  quaternions  can  be  treated  as  simply  the  result  of 
rendering  down  redundant  matrix  information  while  retaining  “everyday”  multiplication.  In  practice,  it’s 
tedious  to  multiply  quaternions  using  i,  j,  k;  instead  we  multiply  them  using  a  non-standard  multiplication. 
But  aside  from  that,  quaternions  still  replace  a  9-element  matrix  by  a  4-element  array,  so  have  advantages 
of  economy  in  numerical  work. 

4That  is,  e 10  =  1  +  i #  +  (i#)2/2!  +  •  ■  ■  =  cos#  +  Isin#,  where  the  last  equality  comes  about  by  noting 

that  the  real  and  imaginary  parts  of  the  complex  series  are  identical  to  the  Taylor  series  of  cosine  and  sine 

respectively.  Note  that  this  sum  defines  e'e ,  because  there  is  otherwise  no  a  priori  meaning  to  a  complex 

exponential;  the  series  expansion  defines  it  in  the  same  way  as  if  i#  were  real.  But  a  small  amount  of 

analysis  (found  in  books  on  complex  variables)  shows  that  this  definition  gives  elS  all  the  behaviour  that 

we  know  and  expect  of  a  real  exponential — behaviour  such  as  eiee1^  =  e1(-s+<^,  and  so  on. 
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2.2  Analytic  Signals  and  Fourier  Analysis 

Regardless  of  whether  a  real  signal  is  written  as  A(t)  sin[cu0i  +  <j>(t)]  or  A(t)  cos[cu0i  +  4>{t)], 
the  corresponding  complex-number  phasor  A(t)e^Uot^'t^  produced  by  (2.12)  is  called  the 
(corresponding)  analytic  signal.  Its  phasor  rotates  counter-clockwise  with  an  instantaneous 
angular  velocity  of  d/d t  (phasor  angle)  =  cu0  +  which  we  have  assumed  is  not  too 

different  from  cuq  because  the  signal  is  narrow  band. 

This  conversion  of  a  real  signal  to  an  analytic  signal  is  related  to,  but  crucially  not 
the  same  as,  the  idea  of  complex-Fourier  analysing  the  real  signal.  It’s  instructive  to 
compare  these  processes  of  (1)  converting  a  real  signal  to  its  analytic  signal  and  (2)  Fourier 
decomposing  the  real  signal  into  complex  components.  Consider  the  real  signal  cos  cu0i  with 
tug  >  0 — hardly  a  useful  signal,  but  a  sufficient  example  for  this  discussion.  Note  that  a 
real  signal  can  always  be  written  using  only  positive  frequencies;  the  very  idea  of  a  negative 
frequency  emerges  only  in  the  following  Fourier  analysis. 

Converting  the  real  signal  to  its  analytic  signal:  The  analytic  signal  formed  from 
the  real  signal  cos  cu0i  is  a  single  phasor  eX0J°L  Both  the  real  signal  and  its  analytic 
signal  are  considered  to  have  the  same  single  frequency  cuq,  which  is  always  positive. 

Complex-Fourier  decomposing  the  real  signal:  A  real  signal  can  always  be  Fourier 
decomposed  into  sines  and  cosines.  This  decomposition  requires  two  sets  of  summa¬ 
tions  and  integrations,  one  for  the  sines  sin  cut  and  one  for  the  cosines  coscut  (the 
zero  subscript  is  omitted  as  cu  is  an  index  of  summation  or  integration),  so  that  all 
terms  in  the  Fourier  analysis  tend  to  appear  twice.  But  sin  cut  and  coscut  can  both  be 
written  as  linear  combinations  of  the  complex  exponentials  eXU)t  and  e_xaJ<.  Of  course, 
we  would  gain  nothing  by  converting  real  sinusoids  to  complex  exponentials  if  the 
result  was  merely  to  replace  a  “sine,  cosine”  pair  with  an  “e1^,  pair.  But  if 

we  allow  the  frequencies  cu  =  2nf  in  the  complex  Fourier  decomposition  to  assume 
negative  values  as  well  as  positive,  this  pair  of  complex  exponentials  collapses  to  just 
one:  exwt,  where  cu  now  ranges  over  all  real  numbers,  negative  as  well  as  positive, 
in  the  summations  and  integrations.  So  the  Fourier  decomposition  writes  cosc u0t 
(with  cu q  >  0)  as  a  sum  of  two  phasors  1/2  elu°t  and  1/2  e~u°ot  rotating  in  opposite 
directions. 

This  use  of  negative  frequencies  in  Fourier  analysis  might  initially  seem  strange, 
but  they  are  simply  a  way  of  halving  the  amount  of  Fourier  notation  needed — with 
the  added  huge  benefit  of  rendering  the  complex  Fourier  transform  and  its  inverse 
almost  identical,  which  aids  the  mathematics.  However,  this  Fourier  language  now 
allocates  the  idea  of  a  single  frequency  cu  not  to  sin  cut  and  coscut,  but  rather  to  elut. 
So  when  we  Fourier-decompose  the  real  signal  coscu0t  into  two  complex  exponen¬ 
tials,  coscu0t  =  1/2  e1^0*  +  1/2  e"1^0/  we  say  that  coscu0t  contains  two  frequencies:  cu0 
and  — cu0-  The  negative  frequency  — cu0  contains  no  more  information  than  the  pos¬ 
itive  frequency  cu 0,  and  is  present  only  to  simplify  the  Fourier  language.  Perhaps 
these  two  frequencies  cu0  and  — cu0  should  be  called  “Fourier  frequencies”  to  reinforce 
the  idea  that  they  result  from  complex-Fourier  decomposing  a  signal  cos  cu0t  (cu0  >  0) 
that  is  clearly  simple  harmonic  motion  at  a  single  frequency  cu0-  More  usual  is  to 
call  cuq  and  — cu0  simply  the  frequencies,  while  the  positive  cu0  is  called  the  repetition 
rate. 
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Table  1:  A  comparison  of  how  the  basic  functions  are  viewed  in  the  analytic- 
signal  picture  versus  the  complex-Fourier  analysis  picture 

Analytic-signal  picture 

Fourier  picture 

eM,t 

sina^t  or  cost o0t 
General  real  signal 

One  frequency  tv0  >  0 
One  frequency  cj0  >  0 
One  instantaneous 
frequency  tv(t)  >  0 

One  frequency  tv0  >  0  or  <  0 
Two  frequencies,  tv0  and  —  w0 
Frequency  spectrum 
with  tv  =  —00  — >•  +00 

Table  1  summarises  the  difference  between  the  analytic  signal  and  Fourier  decomposing 
the  real  signal.  To  reiterate: 

-  When  converting  the  real  signal  cosu}0t  to  its  analytic  signal  both  the  real 

signal  and  the  analytic  signal  are  considered  to  have  the  same  single  frequency  ujq, 
which  is  always  positive. 

-  Fourier  analysis  of  the  real  signal  cos  u0t  works  with  components  elult  with  u  =  uj0 
and  —tv o,  and  these  components  are  each  considered  to  have  a  single  frequency  tv. 
The  real  signal  cos  tv0t  is  then  considered  to  have  two  frequencies,  one  positive  and 
one  negative. 

The  analytic  signal  results  when  a  real  signal  is  converted  to  a  single  phasor,  and  is 
precisely  the  single-phasor  model  of  an  oscillating  wave  that  is  described  in  any  number 
of  introductory  physics  texts.  In  contrast,  complex-Fourier  decomposing  a  real  signal 
produces  phasors  in  pairs  with  frequencies  of  opposite  sign,  whose  utility  lies  in  their 
simplification  of  signal  processing  mathematics. 

Indeed,  converting  cos  co0t  to  its  analytic  signal  e1^0*  is  indistinguishable  from  singling 
out  the  positive-frequency  term  1/2  elulot  from  the  complex  Fourier  pair  and  doubling  its 
amplitude.  Because  of  this,  signal  processing  textbooks  sometimes  describe  creating  the 
analytic  signal  as  “discarding  the  negative  frequency”.  Unfortunately,  this  mixes  complex- 
Fourier  language — which  considers  coscuq t  to  have  two  frequencies  of  opposite  sign — with 
analytic-signal  language,  which  recognises  only  positive  frequencies  and  so  considers  cos  iv0t 
to  have  a  single  frequency.  The  analytic  signal  in  principle  has  nothing  to  do  with  discard¬ 
ing  negative  frequencies,  and  yet  it  can  be  constructed  in  practice  by  discarding  negative 
frequencies  and  doubling  the  weights  of  the  positive  frequencies.  This  is  explained  further 
in  Appendix  A. 3. 

The  signals  emitted  by  most  radars  are  not  pure  sinusoids,  but  they  tend  to  be  suf¬ 
ficiently  narrow  band  that  a  phasor  can  still  meaningfully  be  used  to  portray  them.  In 
this  report  we  assume  that  the  conversion  of  the  received  real  signal  to  an  analytic  signal 
(the  production  of  I/Q  data)  has  already  occurred  in  the  radar  receiver,  so  that  we  deal 
purely  with  the  analytic  signal  at  all  times.  That  is,  the  radar  receiver  can  be  viewed  as 
converting  the  incoming  real  signal  to  a  stream  of  complex  numbers.  We  also  write  the 
signal  just  prior  to  being  emitted  as  a  complex  number,  remembering  that  in  practice  the 
mixing  operation  in  (2.7)  will  then  be  used  by  the  radar  hardware  to  convert  that  complex 
number  to  a  real  signal,  which  is  then  physically  emitted. 
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In  general  of  course,  a  real  signal  s(t)  could  be  converted  to  any  one  of  an  infinite 
number  of  different  complex  signals:  just  add  an  imaginary  part  which  can  be  anything. 
But  only  one  of  these  complex  signals,  the  analytic  signal  sc(t),  is  considered  to  be  special, 
due  to  its  useful  properties  discussed  in  Appendix  A.  That  appendix  describes  the  Hilbert 
transform,  a  procedure  that  constructs  the  analytic  signal  from  a  given  real  signal  after 
the  entire  real  signal  has  been  received.  This  is  in  contrast  to  the  mixing  of  the  real  signal 
with  sin  tag t  and  cos  uj^t  that  is  done  by  a  radar  receiver  while  it’s  receiving  that  signal. 

If  the  real  signal  is  narrow  band  around  a  carrier  frequency  ujq,  the  analytic  signal 
will  contain  a  factor  of  e^0*.  This  mathematical  factoring  out  of  the  zero- information 
carrier  from  the  information-carrying  modulation  renders  some  common  signal  processing 
manipulations  easier  and  tidier.  For  example,  consider  an  emitted  carrier  wave  with 
frequency  u ;0  that  is  bounced  from  a  receding  target  and  returns  “red-shifted”  asw0  —  2  Hz. 
If  we  were  to  use  real  sinusoids  in  our  analysis,  we  would  compare  the  emitted  wave’s 
time- varying  part  sincugf  with  the  received  time- varying  part  sin(u;o  —  2Hz)f.  To  extract 
the  Doppler  frequency  and  its  sign,  we’d  need  to  refer  always  to  the  carrier  w0,  and 
would  note  that  the  received  frequency  was  less  than  the  carrier  by  2  Hz,  so  that  the 
Doppler  frequency  was  —2  Hz,  corresponding  to  a  receding  target.  On  the  other  hand, 
complex  notation  renders  the  extraction  of  the  Doppler  shift  much  more  straightforward: 
we  have  an  “emitted  signal”  e^0*  and  a  “received  signal”  e^w°— 2Hz^,  and  one  stage  of 
our  mathematical  analysis  simply  removes  the  carrier  by  multiplying  the  returned  signal 
by  e-1^0*,  which  we  can  do  in  principle  by  simply  ignoring  e1^4  in  the  mathematics.  The 
result  is  e^~2Hz)4,  which  trivially  yields  an  angular  Doppler  frequency  of  —2  Hz.  So  the 
Doppler  frequency  and  its  correct  sign  emerge  immediately  in  the  complex  approach. 

2.3  Visualising  Transmission  and  Reception  of  a  Signal 

The  above  discussion  makes  the  point  that  when  a  radio  signal  enters  a  radar  receiver, 
the  receiver  produces  a  stream  of  complex  numbers.  It’s  not  the  case  that  we  are  merely 
using  complex  numbers  as  a  redundant  way  to  represent  real  numbers  for  mathematical 
convenience.  Instead,  after  processing,  the  received  signal  has  been  converted  to  a  stream 
of  complex  numbers,  with  no  redundancy  in  its  information  content. 

With  this  in  mind,  we  begin  by  constructing  the  signal  received  by  a  radar  whose 
transmitter  and  receiver  are  co- located  (as  they  always  are  in  this  report).  Write  the 
narrow-band  analytic  signal  emitted  at  time  t  at  the  transmitter/receiver  site  as 

transmitted  signal  sc(t)  =  u(t )  eXUJot ,  (2-13) 

where  u(t)  is  the  complex  form  of  the  I/Q  data:  the  signal’s  modulation.  The  signal  only 
attains  this  far- field  form  at  a  distance  from  the  radar  of  several  wavelengths,  by  which  time 
the  near-field  contribution  has  become  insignificant.  But  that  distance  is  small  enough 
that  we  can  consider  the  signal  to  be  sc(t)  at  the  radar  itself. 

Suppose  that  this  signal  is  bounced  off  a  target  at  a  distance  r,  which  may  be  moving 
with  a  constant  radial  velocity5  v  =  dr/dt.  We  ask:  what  analytic  signal  gc(t)  is  received 
by  the  radar  at  time  f?  (Note  that  the  word  “analytic”  is  understood  as  always  present 
and  will  be  omitted  from  now  on.) 

5Note  that  v  is  positive  for  a  receding  target.  Some  radar  books  use  this  convention,  while  others 
define  v  =  —  dr/dt,  making  v  negative  for  a  receding  target. 
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c  x  time 


Figure  1 :  A  plot  of  c  x  time  versus  radial  distance  for  the  motion  of  both 
radar  and  target,  where  c  is  the  speed  of  light.  The  scenario  is  analysed  in  the 
radar’s  frame,  so  the  radar  is  at  r  =  0  for  all  time.  The  target  is  at  r  =  r0 
at  time  0.  For  a  receding  target,  v  >  0. 


The  scenario  is  shown  in  Figure  1.  With  c  the  speed  of  light,  this  figure  plots  c  x  time 
versus  radial  distance  for  the  motions  of  radar  and  target,  and  incorporates  the  radar  sig¬ 
nal.  Because  the  vacuum  speed  of  light  is  always  measured  by  a  constant-velocity  emitter 
or  receiver  to  have  the  same  value  (with  corrections  for  air  and  a  possibly  accelerated 
frame  being  negligible  here),  we  can,  with  complete  generality,  analyse  the  scenario  in  the 
radar’s  frame,  meaning  we  can  set  the  radar  emitter/receiver  to  be  at  rest  throughout  the 
scenario.  That  is,  only  the  relative  motion  of  radar  and  target  matters,  so  place  the  radar 
at  r  =  0  for  the  duration  of  the  scenario,  with  the  target  at  r  =  r0  at  t  =  0.  To  facilitate 
later  discussion,  suppose  that  t  =  0  starts  a  coherent  processing  interval ,  to  be  studied 
later  in  Section  5. 

One  very  useful  feature  of  this  diagram  is  that  because  A r  =  cAt  for  light,  the  radar 
signals  always  have  a  45°  slope,  which  allows  us  to  make  quick  use  of  simple  geometry  to 
analyse  the  interactions  presented  on  it. 

The  signal  emitted  at  some  general  time  t  is  sc(t)  =  u(t)eluJ°t.  What  is  received  at 
some  general  time  t?  This  received  signal  is  essentially  just  a  lower-amplitude  version 
of  what  was  emitted  at  the  earlier  time  t  —  td,  where  td  is  the  time  delay  from  emission 
to  reception.  The  signal  emitted  at  t  —  td  was  u{t  —  td)  .  After  bouncing  from 

the  target,  this  signal  is  received  almost  unchanged  at  time  t  except  for  the  presence  of 
(a)  some  amplitude  A(rbounce)  that  depends  on  the  target’s  distance  rbounce  at  the  time 
of  bounce,  and  (b)  a  possible  phase  shift  generated  when  the  signal  interacted  with 
the  target,  which  we’ll  absorb  into  A(rbounce)  by  making  that  amplitude  complex.  So  the 
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signal  received  at  time  t  is 


9c(t)  =  ^ (^bounce)  «(*  “  U)  eW°(t  *d).  (2.14) 

The  geometry  of  45°  light  rays  in  Figure  1  indicates  immediately  that 

ctd  2rbounce.  (2.15) 


Thus  (2.14)  becomes 


gc(t)  =  A(rhounce)  u(t  -  td)  e~^o ^houncJc_  (2T6) 

The  target  range  at  the  time  of  bounce  is  rbounce  =  r0  +  v(t  —  td/ 2).  How  does  this  value 
compare  with  the  range  at  the  time  of  emission,  +  v(t  —  td),  and  the  range  at  the  time 
of  reception,  r0  +  vt?  The  difference  in  both  cases  is  vtd/ 2.  What  is  td?  Since  the  target 
travels  distance  a  at  velocity  v  in  time  t  —  td/2,  we  know  that 

ctd/2  =  r0  +  a,  and  a  =  v(t  —  td/2) .  (2-17) 


It  follows  that 

ctd  =  ^r°  +^  ~  2(r0  +  vt) ,  (2.18) 

1  +  v/c 

where  the  last  approximation  holds  because  generally  the  target  speed  |u|  will  be  much 
less  than  c.  Then 

vtd/2  =  v/c  x  ctd/2  ~  v/c(r0  +  vt) .  (2.19) 

Suppose  r0  =  104  m,  v  =  10  m/s,  and  the  current  pulse  is  32  pulses  into  a  coherent  pro¬ 
cessing  interval.  The  time  from  the  start  of  one  pulse  to  the  next  is  100  /xs,  so  that 
i  ~  32  x  100  fi s.  The  range  difference  vtd/2  is  then,  making  use  of  SI  units, 

vtd/2  ~  x  (104  +  10  x  0.0032)  metres  ~  0.5  mm.  (2.20) 

Note  also  that  vt  —  3  cm  here,  so  \vt\  <S  r0,  simplifying  (2.19)  to 


td  ^  2 r0/c. 


(2.21) 


The  difference  vtd/2  between  the  ranges  at  emission,  bounce,  and  reception  is  so  small 
that  we  are  free  to  replace  the  range  at  bounce  in  (2.16)  with  either  the  range  at  emission 
or  reception.  But  it’s  important  to  note  that  comparing  (2.15)  with  (2.21)  doesn’t  allow 
us  to  replace  rbounce  with  the  range  at  the  start  of  the  coherent  processing  interval  r0 
when  Doppler  measurements  are  involved,  because  that  replacement  will  discard  all  infor¬ 
mation  about  the  target  velocity,  which  won’t  allow  later  Doppler  processing  to  be  used. 
Even  so,  the  amplitude  of  the  pulse  is  not  sensitive  to  the  changing  ranges  involved,  so 
can  be  calculated  at  the  start  of  the  coherent  processing  interval,  meaning  for  range  r0. 
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Equation  (2.16)  becomes 

ig—iojQ2/c  x  range  at  emission  'j 
g— ioJ02/c  X  range  at  bounce  l  (any  c]-10Jce  jg  valid) 
x  range  at  reception  I 

!e-tw02/c[r-0  +v(t-td)\ 
e-iu02/c[r0+v(t-td/2)] 
e-iui02/c[r0+vt] 

!e-iiu02  v(t-td)/c  \ 

e-iu02v{t-td/2)/c  l  _  (2>22) 

e~iui02vt/c 

The  first  option  in  (2.22)  is  the  simplest  when  writing  modelling  code,  and  last  option  is 
the  simplest  to  analyse  mathematically: 

Received  signal  when  used  for  modelling: 

gc(t)  ~  A(r0)  u(t  -  2 r0/c)  eiuJot  e~iuJ°2/c  x  range  at  emission 

I _ _ _ I 

time  of  emission 

Received  signal  when  used  for  analysis: 

gc(t)  ~  A(r0)  u(t  -  2 r0/c)  e ^  e~^r0/c  e~iu02vt/c 


(2.23) 


(2.24) 


The  Doppler  Frequency 

Consider  that  a  wave’s  angular  frequency  is  the  rate  of  increase  of  its  phase  at  a  fixed 
position:  uj  =  dphas e/dt.  (Refer  to  (7.6)  ahead  for  more  analysis  of  this  last  expression.) 
So,  from  (2.24),  the  angular  frequency  of  the  returned  wave  measured  at  the  receiver  is 
ur  =  d/dt  (u0t  —  uj02r0/c  —  oj02vt/c).  When  the  radial  target  velocity  v  is  constant  (which 
it  is  to  a  high  approximation  over  the  course  of  the  radar  interaction6),  the  received 
frequency  is 

ur  ~  oj0  —  ujq2v/c  .  (2.25) 

The  (angular)  Doppler  frequency  up  is  defined  to  be  the  increase  in  carrier  frequency  from 
transmission  to  reception: 

ujd  =  ur  —  oj0  ~  —  Uq  2v/c  .  (2.26) 

The  usual  (i.e.  non-angular)  Doppler  frequency  follows  from  (2.26): 

fD  =  ^  =  -fo2v/c.  (2.27) 

This  is  often  written  as  fjj  =  —2v/X0,  with  A0  the  wavelength  corresponding  to  the  carrier 
frequency  /q.  All  of  the  expressions  above  assume  a  narrow-band  modulation  u(t)  that, 
to  a  close  approximation,  doesn’t  appreciably  alter  the  carrier  frequency.  This  is  certainly 
the  case  with  radar  signals. 

6When  v  is  not  constant,  (2.26)  acquires  a  term  —  u>0  2vt/c,  which  will  generally  be  very  small. 
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Calculating  the  Amplitude  A(t) 


We  find  the  amplitude  A  by  asking  how  the  amplitude  of  a  radar  wave  relates  to  its 
transmitted  power.  The  standard  quantity  in  electromagnetic  theory  is  the  Poynting 
vector  S,  which  has  the  direction  of  energy  flow  of  the  wave,  and  whose  magnitude  is  the 
wave’s  areal  power  density,  being  its  power  per  unit  area  normal  to  that  direction: 

S  =  £0c2E  x  B  ,  (2.28) 

where  e0  is  the  permittivity  of  the  vacuum  (we  treat  air  as  approximately  a  vacuum  as 
far  as  radar  propagation  is  concerned),  and  c  is  the  speed  of  light  in  a  vacuum.  The 
electric  and  magnetic  fields  of  a  vacuum  electromagnetic  wave  are  always  perpendicular 
with  strengths  related  by  E  =  cB,  so  it  follows  that  the  magnitude  of  the  Poynting  vector 
is  S  =  £0cE2 .  The  value  of  l/(e0c)  is  approximately  377  U,  the  famous  “impedance  of  the 
vacuum”. 

Far  from  the  source,  the  radiated  electromagnetic  field  can  be  visualised  as  a  sinusoid 
of  amplitude  A,  so  its  electric-held  strength  is  E  =  Asin(cj0t  +  <£>).  We  see  that  the  areal 
power  density  S  =  £0c  E 2  fluctuates  sinusoidally  over  one  period  of  the  carrier  wave;  this  is 
too  rapid  to  be  measured,  but  it’s  sufficient  to  consider  the  average  areal  power  density  (S) 
over  a  period  of  the  wave: 

(S)  =  £0c  {E2)  =  £0c  (A2  sin2(w0t  +  4>))  =  £qcA2/2  .  (2.29) 


Clearly,  a  similar  expression  holds  if  the  wave  is  represented  by  a  complex  exponential 
with  complex  amplitude  A:  just  replace  A2  by  |A|2.  If  u(t)  has  value  0  when  the  radar 
is  not  emitting  a  pulse  and  modulus  1  when  the  radar  is  emitting  a  pulse,  we  can  calcu¬ 
late  A  from  (2.29).  First  consider  one-way  propagation.  Here  subscript  T  denotes  target, 
t  denotes  transmitted,  r  denotes  received  by  the  radar,  and  we  require  ATl  the  value  of  A 
immediately  before  the  wave  strikes  the  target. 


Aj1 


2 


so  that 


—  x  average  areal  power  density  received  by  target  at  r 

£0C 


Pr 


£0c  target  area 


= 

£0C  Gt  A2 


2  PGA2 


4-7T 


£qC  (4-7rr)2  T  Gt A2 


AG, 

2tt£0c  r2  ’ 


Prp  from  one-way 
radar  equation 


At 


1JPtGt  'e1^ 

r  V  27T£qC 


(2.30) 

(2.31) 


for  some  phase  (j>0. 

For  two-way  propagation,  the  value  Ar  of  A  for  the  signal  received  after  being  bounced 
from  the  target  is  not  simply  given  by  (2.31)  with  r  replaced  by  2 r,  because  the  interaction 
with  the  target  complicates  the  picture: 


I AJ2  = 


£0C 

2 


x  average  areal  power  density  received  from  target  at  r 


P, 


x 


£qC  receiver  area 

^PtGta 
£0c  (47r)2  r4  ’ 


£0c 


p. 


47T 

GA2 


2  PtG2a\2  4v r 


£qC  (47t)3  r4  Gt  A2 

i _ i 

Pr  from  two-way 
radar  equation 


(2.32) 
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so  that 


Ar 


1  j2PtGt^ci^ 

47T r2  V  £qc  ’ 


(2.33) 


where  this  is  in  general  different  to  the  0O  in  (2.31).  We  are  interested  in  two-way 
propagation,  so  focus  on  (2.33).  Equations  (2.23)  and  (2.24)  become 


Received  signal  when  used  for  modelling: 


9c(t)  ^ 


1  l2PtGta 


e0c 


A4>o  u(t  —  2r0/c)  eiulot  e~iuj°2^c  x  range  at  emission  (2.34) 


time  of  emission 

Received  signal  when  used  for  analysis: 


9c(t)  ^ 


1  2PtGta 


47rrQ 


£0C 


2^°  u(t  -  2 r0/c)  e^ot  e~iuo2ro/c  e^nt  (2.35) 


As  expected,  the  power  in  the  received  signal  (which  is  proportional  to  its  amplitude 
squared)  is  proportional  to  the  fourth  power  of  target  distance,  the  effective  radiated 
power  PtGt ,  and  the  target  cross  section  a. 

The  determiner  of  range  in  (2.35)  is  the  time  delay  td  =  2 r0/c  from  transmission  to 
reception.  Later  we’ll  discuss  the  correlation  procedure  that  extracts  this  time  delay,  and 
hence  the  range  r0.  Equation  (2.35)  also  includes  a  constant-phase  term  of  and 

of  course  the  unknown  phase  e1^0  introduced  when  the  signal  interacted  with  the  target. 
It  doesn’t  matter  that  the  radar  receiver  cannot  be  aware  of  these  additional  phases, 
because  any  constant-phase  term  in  (2.35)  won’t  affect  the  range-Doppler  calculations  in 
Section  4.2  and  beyond. 


3  Modelling  the  Returned  Signal 

For  the  sake  of  visualisation,  we’ll  assume  the  radar  illuminates  a  ship  that  is  represented 
by  a  large  number  of  known  scatterers.  So  we  require  to  model  the  signal  returned  from 
these  scatterers,  given  some  known  emitted  signal.  If  the  emitted  and  received  signals  are 
represented  by  sampled  I/Q  data,  two  ways  to  model  the  received  signal  are  as  follows. 

The  high-fidelity  way:  We  can  construct  the  received  signal  at  all  sampled  times,  by 
back-tracing  the  rays  that  enter  the  receiver  to  their  sources,  one  at  each  scatterer. 
This  method  works  perfectly  well  but  is  very  slow,  because  it  requires  the  emitted 
signal  to  be  calculated  many  times.  For  example,  if  there  are  1000  scatterers  and  500 
sampling  times  at  which  we  require  the  received  signal,  then  at  each  sampling  time, 
we  must  trace  the  signal  via  1000  scatterers  back  to  the  emitter,  thus  calculating  1000 
values  of  the  emitted  signal  for  each  of  those  500  times,  making  500,000  calculations 
of  the  emitted  signal.  This  approach  is  prohibitively  CPU-expensive. 

The  slightly  lower- fidelity  way:  This  approach  makes  use  of  linear  time-shift-invariant 
signal  processing  theory.  It  only  requires  the  returned  signal  to  be  calculated  corre¬ 
sponding  to  a  ping  signal  emitted  at  time  zero — meaning  that  a  total  of  only  1000 


UNCLASSIFIED 


13 


DSTO-TN-1386 


UNCLASSIFIED 


m 


LTSI  system 


h[t]  known 

- >- 


m 


LTSI  system 


g[t]  =  ? 

- 


Figure  2:  A  “black  box”  that  illustrates  the  approach  of  calculating  the 
response  of  an  LTSI  system.  Top:  given  an  impulse  <5[£],  we  know  the  re¬ 
sponse  h[t\.  Bottom:  given  a  signal  f[t\,  what  is  the  response  g[t\? 


calculations  are  needed  for  the  above  1000  scatterers.  A  further  processing  gain  can 
also  be  achieved,  as  some  of  the  calculations  depend  only  on  the  emitted  signal  and 
so  can  be  made  off  line.  We’ll  use  this  approach  in  this  report  and  will  describe  it 
in  detail  now. 


Linear  time-shift-invariant  theory  calculates  the  output  of  a  “linear  time-shift-invariant 
(LTSI)  system”  given  some  input,  subject  to  two  assumptions: 


-  the  system  is  linear :  if  we  e.g.  triple  the  strength  of  the  input  signal,  the  strength 
of  the  output  signal  will  be  tripled;  and  if  two  signals  are  input  simultaneously,  the 
output  will  be  the  sum  of  the  individual  outputs  from  those  signals; 

-  the  system  is  time-shift  invariant:  whether  a  signal  is  input  today  or  tomorrow 
doesn’t  affect  the  output  signal.  Clearly  this  only  applies  to  a  moving  ship  over 
short  time  intervals. 


The  radar  interaction  is  certainly  linear.  As  for  time-shift  invariance,  how  small  must  the 
relevant  time  interval  be  for  the  invariance  to  apply?  We  saw  in  Figure  1  and  (2.19)  that 
an  emitted  signal  will  return  after  a  time  interval  t ^  ~  2 (tq  +  vt)/c.  Only  when  vt  <C  r0 
will  there  be  no  “stretching”  of  flight  times  of  successive  signals.  This  inequality  holds  for 
the  typical  example  considered  in  Section  2.3  so,  for  our  scenarios,  the  time-shift  invariance 
holds  well. 

The  central  result  of  LTSI  theory  is  that  the  signal  output  from  an  LTSI  system  equals 
the  convolution  of  the  signal  input  with  the  system’s  impulse  response.  The  system’s 
impulse  response  is  its  output  when  the  input  is  an  impulse:  a  “spike”  in  the  continuous 
case,  or  a  “1”  in  the  discrete  case.  We  prove  this  central  result  while  referring  to  Figure  2. 
Consider  a  signal  that  has  been  sampled  at  times  t  to  give  numbers  f[t\  (the  brackets 
denote  the  discreteness  of  /).  It  enters  the  LTSI  “black  box”  and  is  transformed  into  an 
output  g,  sampled  as  g[t\,  where  we  denote  the  transformation  as  g  =  L(f).  Given  f[t\, 
what  is  g[t]7  We  know  the  system’s  response  h[t]  when  a  single  “1”  is  input  to  it  at 
time  zero.  This  signal  that  is  “1”  at  time  zero  and  “0”  at  all  other  times  is  denoted  <5[£]; 
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that  is,  <5[0]  =  1  and  5[t  /  0]  =  0.  So  we  know  that  h[t]  =  L(S[t]).  Now: 


»M  =  u/[(])  =  r(£ /[(>[*-*']) 

t> 


lincarity  £)/[*']  £(<5[t-t']) 


time-shift 

invariance 


t> 


=  (f*h)  [t] , 


where  f  *  h  is  the  convolution  of  /  and  h: 

(f  *  h)  [t]  =  Y,mh[t-t']. 

t' 


(3.1) 


(3.2) 


The  key  equation  of  LTSI  theory  is  g  =  /  *  h,  or 


(3.3) 


Appendix  B  has  a  tutorial  on  convolution;  I  think  that  practice  with  convolving  two  arrays 
from  first  principles  using  pen  and  paper  is  invaluable  to  understanding  the  procedure. 


output  signal  =  input  signal  *  impulse  response . 


3.1  Calculating  the  Ship’s  Impulse  Response 

To  calculate  the  elements  of  the  ship’s  impulse  response  h[t],  realise  that  they  are  the 
returns  of  an  emitted  pulse  that  is  a  single  “1”.  That  is,  emit  a  single  “1”  at  time  zero 
and  place  the  returns  from  all  the  scatterers  (in  temporal  order  of  course)  into  an  array  h, 
which  becomes  the  impulse  response  for  use  in  later  calculations.  Of  course,  we  don’t 
assume  that  we  know  where  the  ship  is,  and  so  we  should  calculate  the  impulse  response 
of  a  region  of  interest:  everything  that  lies  between  some  range  RneaT  to  some  range  i?far. 
This  could  include  sea  clutter,  but  in  the  following  discussion  I’ll  assume  without  loss  of 
generality  that  only  the  ship  gives  a  return. 

The  impulse  response  h  of  the  “system”  (ship  plus  environment)  is  supposed  to  be 
everything  returned  from  the  moment  that  we  emit  the  “1”  impulse.  That  is,  we  emit 
an  impulse  and  then  immediately  listen  to  hear  what  the  world  throws  back  at  us.  We 
can  begin  to  create  this  impulse  response  by  initialising  an  array  of  zeroes,  some  of  which 
will  eventually  be  replaced  by  returns.  How  many  zeroes  are  in  this  array?  As  shown  in 
Figure  3  on  the  following  page,  the  array  is  initialised  with  zeroes  that  are  place  holders 
for  samples  taken  at  time  intervals  of  ts.  In  principle  we  then  calculate  the  return  of  each 
scatterer  in  the  world,  and  add  that  return  to  the  nearest  element  in  time  of  the  impulse 
response  array  h.  But  in  practice  we  need  only  calculate  the  returns  of  the  scatterers 
in  the  region  of  interest.  This  means  that  the  red  zeroes  in  the  figure  don’t  interest  us, 
and  for  efficient  memory  management  we  can  ignore  them,  and  account  for  their  absence 
simply  by  adding  /?,near  to  all  ranges  that  are  calculated  from  the  impulse  response  that 
was  initialised  with  the  black  zeroes  in  Figure  3. 

So  given  a  sampling  interval  of  ts,  the  impulse  response  is  initialised  to  an  array  of 
M  zeroes,  such  that  (M  —  l)fs  is  the  time  interval  between  the  return  from  i?near  to  the 
return  from  i?far: 

(M  -  1  )ts  ~  2 Riar/c  -  2 Rneai/c .  (3.4) 
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Figure  3:  The  idea  behind  calculating  the  ship’s  impulse  response:  create 
an  array  of  zeroes,  then  increment  its  relevant  elements  by  the  returns  from 
a  “1  ”  sent  out  that  bounces  from  each  scatterer. 


In  practice  we  must  ensure  that  M  is  a  whole  number: 


M  —  1  =  round 


2(-Rfar  —  -^near  ) 


ct. 


(3.5) 


Now  bounce  a  signal  “1”  off  each  scatterer  and  find  the  time  when  it  returns,  adding  the 
returned  signal  to  the  element  of  h  whose  index  is  the  nearest  one  that  corresponds  to 
that  time  of  return.  For  scatterer  n  at  range  rn  the  nearest  index  will  be  kn,  where 


—  2rn/c  2I?near/  c , 


again  with  some  rounding: 


kn  —  1  =  round 


cts 


(3.6) 


(3.7) 


The  return  of  the  impulse  “1”  is  the  return  of  a  signal  with  modulation  1  emitted  at 
time  zero.  Now  use  (2.34)  to  write 


modulation  of  return 


4vr rl  V  e0c 


icj02/c  x  range  at  emission 


(3.8) 


Calculate  this  quantity  for  each  scatterer;  e.g.,  calculate  it  for  scatterer  n  and  add  this  to 
element  kn  of  the  array  h,  found  from  (3.7).  This  modulation  must  be  added  to  element  kn 
because  it  might  well  be  that  two  closely  spaced  scatterers  (in  range)  both  contribute  to 
this  element. 
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4  The  Correlation  Procedure 

When  our  radar  emits  a  signal  that  reflects  from  a  stationary  target  and  returns  to  us  in  a 
degraded  form,  we  can  correlate  the  emitted  and  received  signals  to  calculate  the  round- 
trip  time.  Each  signal  is  sampled  and  digitised.  Write  the  sampled  emitted  signal  (carrier 
removed)  as  an  array  E  of  complex  numbers,  of  length  e.  Likewise,  the  received  signal 
(sampled  at  the  same  rate  and  with  carrier  removed)  is  an  array  R  of  length  r.  Further, 
denote  the  complex  conjugate  of  a  number  z  by  z* ,  and  introduce  the  idea  of  the  complex 
conjugate  of  a  sequence  and  its  reverse.  For  example,  a  3-sample  array  E  =  [Ey  E2  E3] 
has 

E*  =  [E\  El  El]  ,  and  E]  =  [E|  El  E{]  .  (4.1) 

Write  the  correlation  of  E  with  R  as  a  new  sequence  E  *  R. 

Calculating  E  *  R  is  visualised  in  (4.2)  below.  Write  the  elements  of  E*  as  a  length-e 
movable  row  above  a  stationary  row  of  the  r  elements  of  R,  such  that  the  last-received 
(rightmost)  element  of  E*  lies  above  the  first  element  of  R,  and  then  multiply  these  two 
elements;  this  product  E^Ry  is  the  first  element  of  E  *  R.  To  calculate  the  next  element 
of  E  -k  R,  shift  E*  one  element  to  the  right  and  again  multiply  each  element  of  E*  with 
the  element  of  R  (if  any)  immediately  below  it.  There  will  now  be  two  such  products,  and 
they  are  added  to  give  the  second  element  of  E  *  R.  Continue  this  procedure  until  there 
are  no  more  products  to  be  calculated.  In  other  words,  after  each  shift  of  E*  one  element 
to  the  right,  we  form  the  dot  product  of  the  overlapping  subsequences: 

E*  =  [E{  El  El  ...  E*e\  — moving  to  right 

R  =  [  R\  R2  R3  ■  ■  •  Rr  \ 

received  signal 


E  *  R  =  [  E%RU  El_\R\  +  E*i?2,  El_2Rx  +  E%_XR2  +  E*eR3,  ...,  ElRr]  (4.2) 

[The  correlation  process  is  sometimes  defined  using  E  and  R* ,  which  differs  only  by  a 
complex  conjugate  from  (4.2)  and  so  is  just  as  valid  to  use.]  This  procedure  will  always 
produce  a  peak  in  the  absolute  value  of  the  correlation  when  a  sequence,  on  being  correlated 
with  itself,  matches  up  with  itself — which  is  the  whole  point  of  correlating  two  sequences: 
to  search  for  one  within  the  other.  (The  proof  is  omitted  but  straightforward,  and  uses 
the  Cauchy-Schwarz  theorem  for  complex  sequences.)  This  peak  would  not  be  guaranteed 
to  occur — and  in  practice  will  seldom  occur — if  the  complex  conjugate  were  not  used  in 
the  correlation  procedure. 

Correlation  turns  out  to  be  almost  the  same  process  as  convolution  (not  to  be 
confused  with  a  complex  conjugate  here,  which  uses  this  symbol  as  a  superscript): 

A*JB  =  At*5,  (4.3) 

where  “f”  was  defined  in  (4.1).  This  is  essential  to  note,  because  convolution  can  be 
carried  out  efficiently  using  a  fast  Fourier  transform,  implying  that  correlation  can  be 
implemented  in  the  same  efficient  way.  Of  great  use  is  the  fact  that  convolution  is  com¬ 
mutative:  A  *  B  =  B  *  A,  and  associative:  (A  *  B)  *  C  =  A  *  ( B  *  C ).  Correlation  is 
neither  of  these;  specifically, 

A*B  =  {B*A)\  (A*B)*C  =  At*(B*C).  (4.4) 

Appendix  B  gives  some  further  comments  on  how  correlation  is  related  to  convolution. 
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4.1  Relating  Correlation  to  Distance  and  Velocity 

Each  of  the  above  “correlation  steps”  (shifting  E*  to  the  right  by  one  sample)  corresponds 
to  testing  for  a  new  posited  time  interval  between  emission  and  reception;  this  time  interval 
corresponds  to  a  posited  range  to  the  target.  If  we  plot  the  absolute  value  of  the  correlation 
versus  this  range,  then  the  range  value  at  which  the  graph  peaks  is  where  we  presume  the 
target  to  be. 

The  target  may  well  be  moving,  in  which  case  we  cannot  expect  a  good  correlation  of 
emitted  and  returned  signals.  Instead,  we  must  allow  that  R  might  be  Doppler  shifted 
from  E.  This  might  be  arranged  by  performing  the  above  correlation  process  many  times; 
each  time,  we  subtract  some  amount  of  posited  Doppler  frequency  from  R.  What  results 
is  many  plots  of  correlation  output  versus  range:  one  plot  for  each  of  the  posited  Doppler 
frequencies.  These  plots  can  be  stacked  to  form  a  surface  plot  of  correlation  output  (or 
rather  its  absolute  value,  since  it’s  generally  complex)  versus  range  and  Doppler  frequency. 
The  properties  of  this  graph,  the  ambiguity  function  for  the  radar  signal  used,  have  been 
studied  extensively,  and  it  forms  the  standard  approach  to  determining  how  well  any  newly 
designed  signal  might  perform  in  practice.  We  might  hope  to  find  a  sharp  spike  in  the 
function  at  the  correct  values  of  range  and  Doppler  frequency.  In  fact,  one  of  the  various 
theorems  of  the  function’s  properties  discussed  in  [2]  states  that  it’s  impossible  to  craft  a 
signal  that  will  produce  such  a  spike.  Appendix  G  shows  how  the  ambiguity  function  is 
created  in  practice. 

The  ambiguity  function  is  built  by  stacking  many  correlation-range  plots,  each  for  a 
particular  value  of  Doppler  frequency.  Historically,  radar  receivers  have  not  calculated 
target  range  and  velocity  in  this  way,  perhaps  because  of  the  large  amount  of  computer 
processing  required.  Instead  they  build  a  range-Doppler  plot  by  stacking  many  frequency 
spectra  of  correlation  data  that  have  been  obtained  pulse  by  pulse  from  a  coherent  pulse 
train ,  with  each  spectrum  calculated  for  a  particular  value  of  range.  The  details  of  this 
procedure  are  given  in  Section  5. 

Summarising,  the  correlation  of  emitted  and  received  signals  is  traditionally  done  in 
either  of  two  ways: 

The  ambiguity  function:  For  each  value  of  Doppler  frequency  in  an  appropriate  do¬ 
main,  plot  the  correlation  of  part  or  all  of  the  signal  sent  out  versus  range,  and 
stack  the  graphs  across  frequency.  This  approach  is  used  off-line  for  analysing  radar 
waveforms. 

The  range— Doppler  plot:  For  each  value  of  range  in  an  appropriate  domain,  plot  the 
frequency  spectrum  of  the  correlations  produced  from  a  series  of  pulses  in  a  coherent 
pulse  train,  and  stack  the  graphs  across  range.  Radar  hardware  does  this. 

The  process  of  correlating  the  returned  signal  with  the  emitted  signal  is  known  as 
applying  a  matched  filter,  whose  name  perhaps  derives  from  something  like  the  following 
argument.  The  primary  process  here  is  that  of  correlating  the  emitted  and  returned  sig¬ 
nals  to  locate  the  emitted  signal  within  the  returned  signal.  In  principle,  this  correlation 
could  be  done  in  an  infinite  number  of  different  ways;  the  question  is,  which  of  these  ways 
can  best  detect  a  signal  in  the  presence  of  noise?  Determining  this  correlation  procedure 
has  traditionally  followed  a  linear  time-shift-invariant  approach  (mentioned  at  the  start 
of  Section  2),  perhaps  because  that  approach  was  historically  the  easiest  method  to  build 
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into  analogue  hardware.  If  the  noise  in  the  returned  signal  is  assumed  to  be  “white” 
(meaning  temporally  uncorrelated),  then  it’s  not  difficult  to  show  [3]  that  the  best  way  to 
correlate  emitted  and  returned  signals  using  a  linear  time-shift-invariant  approach  is  via 
a  convolution.  Convolution  turns  out  to  be  a  type  of  moving  mean ,  and  a  moving  mean 
removes  (“filters”)  signal  components  of  certain  frequencies.  So  in  this  sense  the  corre¬ 
lation  process  can  be  viewed  as  a  type  of  filtering  process,  whose  filter  stencil  “matches” 
the  signal.  Although  the  standard  correlation  procedure  is  often  simply  introduced  as  a 
convolution,  it’s  useful  to  remember  that  convolution  is  just  one  choice  of  how  correlation 
might  be  implemented.  But  convolution  is  the  standard  choice,  because  it  turns  out  to 
maximise  the  ratio  of  what  results  from  filtering  the  signal,  compared  to  what  results  from 
filtering  the  ever-present  noise. 

4.2  Calculating  the  Target  Range 

As  on  page  17,  sample  the  modulation  of  the  emitted  signal  into  a  length-e  array  E  of 
complex  numbers.  Likewise,  sample  the  modulation  of  the  received  signal  at  the  same  rate 
into  the  complex  array  R  of  length  r.  The  sampling  time  interval  is  ts.  We  will  correlate 
E  with  R ,  using  (4.3)  to  produce  an  array  of  correlation  numbers  via  a  convolution. 
Calculating  the  target  range  requires  knowledge  only  of  the  absolute  value  of  each  of  these 
numbers,  so  call  the  array  of  these  absolute  values  C: 

C  =  \E*R\  =  \EUR\,  (4.5) 

which  is  accomplished  by  the  Matlab  command 
C  =  abs(  conv(conj (f liplr (E) ) ,  R)  ); 

The  array  C  has  length  e  +  r  —  1,  as  can  be  shown  by  examining  what  happens  when  E 
(length  e)  is  correlated  with  R  (length  r),  with  reference  to  (4.6)  below.  First,  there  are 
r  correlation  operations  (producing  the  first  r  elements  of  C)  as  the  right-most  element 
of  E*  starts  at  the  left-most  element  of  R  and  shifts  element-by-element  to  eventually 
“hit”  the  right-most  element  of  R.  As  E*  continues  moving  to  the  right,  e  —  1  more 
correlation  operations  occur  until  E*  eventually  “falls  off  the  right-hand  end”  of  R. 

One  of  the  elements  of  C  will  be  a  maximum,  defining  the  best  correlation  of  E  and  R. 
What  time  delay  between  emitted  and  received  signals  does  this  or  any  other  element  of  C 
correspond  to?  For  the  sake  of  the  numerical  argument  below,  suppose  that  the  sampled 
received  signal  contains  the  emitted  signal  somewhere  within  it.  That  is,  suppose  that  R 
is  composed  of  n  noise  numbers  (small  numbers  that  won’t  correlate  well  with  E,  which 
we  replace  by  zeroes  for  the  purpose  of  this  discussion)  followed  by  E  (which  might  have 
some  added  noise,  but  that’s  immaterial  to  this  discussion  so  is  omitted),  followed  by  more 
noise  numbers: 

E*  =  [E{  El  El  ...  E*\  — >  moving  to  right 

R  =  [  0  0  0  ...  0  E1  E2  E3  ...  Ee  0  0  0  ...  0  ] .  (4.6) 

n  noise  numbers  emitted  signal  more  noise  numbers 

The  peak  correlation  occurs  when  E*  has  slid  in  from  the  left  and  now  is  exactly  alongside 
(above)  the  copy  of  E  that  is  inside  R,  element  for  element.  The  number  of  correlation 
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events  that  this  corresponds  to  is  the  sum  of  two  numbers:  n  correlations  with  the  noise 
numbers,  and  then  e  correlations  to  line  up  alongside  E.  So  the  peak  occurs  at  element 
number  n  +  e  of  the  correlation  array  C  (whose  elements  are  numbered  1  to  e  +  r  —  1). 
But  we  know  that  the  time  interval  td  between  emission  and  reception  of  the  signal  is 
composed  of 7  n  lots  of  ts ,  because  there  are  n  —  1  inter-element  occurrences  of  ts  in  the 
noise  signal,  plus  one  extra  occurrence  of  ts  to  reach  element  Ex : 

o  4^  o  4^  o...o  ^  ...  (4.7) 

i _ 1 1 _ i 

n  noise  numbers,  so  1  lot  of  ts 
n  —  1  lots  of  ts 

i _ i 

td  nts 

Because  the  index  of  the  element  of  C  that  corresponds  to  maximum  correlation  is  n  +  e, 
we  conclude  that 


td  =  nts  =  ( [index  of  peak  correlation  of  C  ]  —  e)  ts  .  (4-8) 

(This  equation  refers  to  the  index  of  the  element  of  C  that  gives  the  peak  correlation,  not 
the  element  itself.)  The  time  interval  td  corresponds  to  a  target  range  of  ct^/2,  so  the 
range  inferred  from  the  above  correlation  is 


(4.9) 


As  a  check,  (4.9)  says  that  if  the  correlation  elements  of  C  peak  at  index  e,  the  target 
range  must  be  zero.  This  makes  sense  because  it  corresponds  to  the  absence  of  the  first  n 
noise  numbers  in  (4.6):  the  signal  is  being  returned  as  soon  as  it’s  emitted,  corresponding 
to  a  target  at  zero  distance  from  the  radar. 

Equation  (4.9)  computes  the  target  range  from  the  peak  correlation  index  of  C;  but 
we  can  just  as  well  apply  it  to  every  index  of  C  to  give  the  computed  range  had  that 
index  denoted  the  maximum  element  of  C .  (In  particular,  increasing  the  correlation  index 
by  one  results  in  a  computed  range  increase  of  cts/2.)  This  idea  allows  us  to  rescale  the 
time-delay  axis  to  range,  which  is  a  more  useful  parameter  than  time  delay  for  analysing 
a  real  radar  return. 

Incorporating  the  Target’s  Impulse  Response  h 

Equation  (3.3)  says  that  R  =  E  *  h.  The  correlation  of  emitted  and  received  signals  is 
then 

E*R  =  E]*R  =  E]*(E*h)  =  (E]*E)*h.  (4.10) 

This  is  the  central  equation  that  we’ll  use  to  correlate  emitted  and  returned  signals: 


(4.11) 


7This  calculation  is  “coarse”  by  an  amount  set  by  t3  since  we  are  treating  n  =  td/ts  as  a  whole  number, 
even  though  it  will  not  be  in  practice.  But  when  ts  is  small,  the  range  error  in  the  above  analysis  is 
negligible. 


Correlation  =  E  *  R  =  (E^  *  E)  *  h  . 


ct 

target  range  =  ( [index  of  peak  correlation  of  C }  —  e) 
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The  array  of  emitted  I/Q  data  E  doesn’t  change  from  pulse  to  pulse,  so  E^  *  E  can  be 
pre-computed,  and  then  convolved  at  run  time  with  the  ship’s  impulse  response  h,  which 
is  updated  at  the  start  of  each  new  pulse.  The  absolute  values  of  the  result  form  the 
array  C  in  (4.5)  that  will  return  the  ship’s  distance  via  (4.9).  Also,  the  correlation  arrays 
Ei*  Ri,  E2*R2,  ■  ■  •  from  pulses  1,  2, . . .  form  successive  rows  of  the  matrix  M  discussed 
ahead  on  page  22,  whose  columns  turn  out  to  hold  Doppler  information. 

Note  that  in  general,  prepending  or  appending  zeroes  to  one  or  both  of  any  arrays 
A  and  B  has  the  effect  of  prepending  or  appending  those  zeroes  to  A*  B.  This  fact  is 
worth  remembering  whenever  an  array  begins  with  a  set  of  zeroes  that  we  wish  to  exclude 
for  reasons  of  efficient  calculation.  This  idea  was  actually  implicit  in  the  discussion  on 
page  15  when  we  spoke  of  adding  the  offset  i?near  to  all  ranges  calculated  using  an  h  that 
had  no  unnecessary  leading  zeroes. 

So  much  for  extracting  range  information  from  the  returned  signal.  We’ll  now  show 
how  to  extract  Doppler  information  from  the  correlations. 


5  Creating  a  Range— Doppler  Plot 

A  radar  receiver  generates  a  range-Doppler  plot  by  emitting  a  set  of  N  coherent  pulses, 
then  correlating  each  pulse  with  that  pulse’s  return.  (It  doesn’t  correlate  the  entire  pulse 
train  it  emits  with  the  entire  train  returned;  cf.  the  ambiguity  function  in  Appendix  G.) 
We’ll  see  shortly  that  the  correlations  over  successive  pulses  contain  Doppler  information, 
which  can  be  extracted  for  each  range  cell  from  a  discrete  Fourier  transform  (DFT)  of  the 
correlations  in  that  range  cell. 

For  simplicity,  suppose  we  have  a  single  point  target,  and  consider  the  first  pulse  Ei 
emitted  in  the  coherent  pulse  train.  Its  modulation  u(t)  is  sampled  at  intervals  of  ts  and 
digitised  into  a  sequence  of,  say,  4  complex  numbers: 

Ei  =  [ui  u2  u3  u4]  .  (5.1) 

The  return  Ri  of  the  first  pulse  is  received  beginning  at  some  time  t0.  Referring  to  (2.35) 
for  one  scatterer,  write  i?1  as  some  modulation  times  the  sampled  Doppler  shift: 

Ri  =  [eiulDt °ua,  e^D(t°+ts)u6,  etaJfl(to+2ts)  uc ,  e^c(A+3L) 

=  e^Dt°  [ua  ,  e5"®1*  ub ,  Uc  ?  Ud]  ,  (5.2) 

where  ua,ub,uc,ud  are  values  of  the  signal  emitted  at  earlier  times  that  incorporate  the 
other  phase  shifts  in  (2.35).  For  simplicity  in  what  follows,  write  the  constants  in  these 
expressions  as 

a  =  ei  ^Df0)  fj  =  eWA.,,  (5.3) 

so  that 

E-i  =  [ui  u2  u3  u4]  ,  Rx  =  a[ua  /3ub  f32uc  (33ud]  .  (5.4) 

We  might  expect  the  next  pulse  E2  emitted  in  the  coherent  pulse  train  to  be  identical 
to  Ei,  with  the  target  still  returning  the  modulation  ua,  ub ,  uc,  ud  since  it  moves  negligibly 
in  the  pulse  repetition  interval  T  from  one  pulse  to  the  next.  Or  there  might  be  some 
known  phase  shift  from  one  emitted  pulse  to  the  next,  as  well  as  the  fact  that  the  reference 
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oscillator  used  to  remove  the  carrier  signal  is  not  completely  stable  and  introduces  some 
arbitrary  phase  to  each  emitted  pulse.  So  call  the  total  phase  shift  cpn  for  the  nth  pulse, 
with  c()i  =  0.  This  phase  is  preserved  in  the  returned  pulse: 

E2  =  [u4  u2  u3  u4]  e102 , 

R2  =  [e^n^o+U  ua}  e^n(A+T+L)  ub  ,  e^n(t0+T+2ts)  ^ 

=  ae'“cTK  Pub  fj2'uc  /33Wd]e^2. 

Similarly, 

E3  =  [Ul  u2  u3  u4]  e1^3, 

R3  =  [e™n(t0+2T)  Ua  ,  eiwD(t0+2T+ts)  Ub  ^  eiwD(t0+2T+2ts)  ^ 

=  ae'UJr,2T[ua  (3ub  /32uc  p3ud]el(p3, 

and  in  general,  for  the  nth  pulse  in  the  coherent  train, 

En  =  [«i  u2  u3  u4]  e^"  ,  Rn  =  aex“D(n~1)T[ua  pub  f32uc  f33ud]  e10-.  (5.7) 

Now  correlate  En  with  Rn,  noting  that  the  conjugation  used  in  the  correlation  process 
cancels  the  unknown  phase  </>n: 

En  -k  Rn  *  Rn 

=  \u%  u3  u2  ul]-e^^*  aeUilD<'n~1^T[ua  (3ub  (32uc  /33ud]je^ 

=  aeiUD<'n~1)T[ulua  ,  ulj3ub  +  u*3ua,  ...,  u\p3ud]  .  (5.8) 

Now  make  this  sequence  the  nth  row  of  a  matrix  M :  this  matrix  will  have  N  rows  because 
each  row  contains  one  of  the  N  correlations  performed  during  the  coherent  processing 
interval  (CPI)  comprised  of  N  pulses.  The  process  of  correlating  each  of  the  N  pulses 
with  itself  is  illustrated  in  Figure  4.  We  correlate  the  first  pulse  sent  out  with  its  return, 
and  set  the  resulting  sequence  of  numbers  to  be  the  first  row  of  the  matrix  M .  The 
numbers  in  this  row  are  the  correlations  represented  by  the  blue  dots  in  the  top  peaked 
“curve”  of  the  figure.  The  peak  corresponds  to  the  best  correlation,  which  gives  the  range 
to  the  target  using  the  procedure  of  (4.9).  But  we  won’t  actually  use  these  correlation 
numbers  directly  to  calculate  this  range,  because  after  having  done  the  correlations,  it 
turns  out  that  we’ll  Fourier-transform  the  matrix  M.  However,  the  resulting  matrix  will 
preserve  the  placement  of  the  correlation  peaks. 

After  having  correlated  the  first  emitted  pulse  with  its  return,  we  do  the  same  with  the 
second  pulse  emitted  and  its  return,  setting  the  resulting  sequence  of  numbers  to  be  the 
second  row  of  the  matrix  M.  This  row’s  numbers  are  the  blue  dots  in  the  second  peaked 
curve  of  Figure  4.  (The  peak  will  lie  at  approximately  the  same  range  as  the  peak  in  the 
top  row,  because  the  target  has  barely  moved  in  time  T .)  We  do  this  for  all  N  pulses — all 
of  the  blue  dots  in  the  figure — to  build  the  entire  correlation  matrix  M,  row  by  row. 

With  M  constructed,  focus  on  any  of  its  columns.  Inspection  of  (5.8)  shows  that  this 


gia;D(t0+2T+3is)^j  gi</>3 

(5.6) 


gicjD(io+^"’+3L),u^J  gi</>2 

(5.5) 
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Figure  4:  The  general  scheme  followed  to  build  a  range-Doppler  plot  via 
range  slices 


column  will  be 


11 


giuD2T 


X 


some  constant  that  changes 
from  column  to  column 


^iu>D(N  —  1)T 


(5.9) 


The  sequence  of  N  numbers  in  this  column  equals  an  irrelevant  constant  times  a  sequence 
of  N  numbers  that  can  be  treated  as  samples  of  at  times  t  =  0,  T,  2 T, ...  So  if  we 
calculate  the  frequency  spectrum  of  the  entire  column,  then — given  that  we’re  using  a 
point  target  in  this  discussion — we  will  (ideally)  find  a  single  angular  frequency  ujjj  to  be 
present.  This  corresponds  to  the  target  recession  velocity  v  via  (2.26)  or  (2.27): 


target  recession  velocity  v 


~cwD 

2w0 


- cfD 
2/o  ' 


(5.10) 


The  frequency  spectrum  of  the  column  in  (5.9)  is  found  using  a  DFT,  which  produces 
an  array  of  the  amounts  of  each  frequency  present.  These  frequencies  are  equally  spaced, 
running  from  approximately  minus  the  Nyquist  frequency  to  the  Nyquist  frequency  for  the 
sampled  set;  see  the  discussion  in  Section  6  for  details  of  the  exact  bounds  of  the  spectrum, 
specifically  (6.4)  and  (6.5).  The  Nyquist  frequency  is  1/(2T),  or  half  the  pulse  repetition 
frequency.  This  array  of  frequencies  can  be  converted  to  radial  velocities  using  (5.10). 
That  is,  each  of  these  frequencies,  if  it  were  the  actual  Doppler  frequency  of  the  target, 
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would  correspond  to  a  velocity  given  by  (5.10),  in  which  case  we  need  only  multiply  the 
entire  array  of  frequencies  by  — c/(2/0)  to  convert  it  from  frequency  to  radial  velocity. 
Ideally,  a  moving  point  target  might  be  expected  to  produce  a  single  sharp  peak  at  its 
Doppler  frequency,  but  in  practice  the  sampling  process  spreads  this  peak  out  over  many 
frequencies. 

We  now  repeat  the  above  DFT  analysis  for  each  column  of  M.  Only  the  common 
factor  outside  the  brackets  in  (5.9)  changes  from  column  to  column  but,  being  constant  in 
time,  it’s  not  detected  by  the  DFT  and  so  doesn’t  affect  the  spectrum  produced.  Hence 
each  column  yields  a  frequency  spectrum,  and  this  column  (with  its  frequency  spectrum) 
corresponds  to  a  particular  range  by  way  of  the  discussion  following  (4.9)  on  page  20.  That 
discussion  is  worth  repeating  here:  if  we  apply  (4.9)  to  every  correlation  index — meaning 
each  column  number  of  M — then  we’ll  obtain  the  range  corresponding  to  that  column 
number,  and  this  converts  the  time  axis  to  range.  So  the  frequency  spectrum  for  each 
column  of  M  becomes  a  slice  of  a  surface  drawn  over  range  and  velocity  axes,  and  this 
surface  is  the  range-Doppler  plot. 

The  same  procedure  applies  to  a  target  composed  of  multiple  point  scatterers.  The 
signals  returned  from  multiple  point  scatterers  simply  add,  and  the  linearity  of  the  convo¬ 
lution  procedure  ensures  that  each  velocity  present  gives  rise  to  its  own  Doppler  peak. 


6  Comments  on  Applying  a  Discrete  Fourier 

Transform 


Different  schemes  for  ordering  the  elements  of  a  DFT  exist.  For  this  report  we  define  a 
length-IV  sequence  {x0, . . . ,xN_i}  to  have  a  length-iV  DFT  of  ,  -X”-[jv/2]+i,  ■  ■  ■  } 

(where,  for  positive  a,  the  notation  [a]  denotes  the  greatest  integer  ^  a).  We  use  the 
convention  of  reference  [4]  with  the  choice  a  =  1  there: 


i 

Xn=^J2  e-i2”kn/N,  n  =  —  [N/2] ,  —  [N /2]  +  1, 
iV  fc  =  o 


(N  values).  (6.1) 


The  inverse  transform  too  is  a  sum  over  N  numbers: 


xk  = 


y  ptenkn/N 


n  =  -[N/  2] 


O^k^N-l. 


(6.2) 


In  practice  for  any  data  set  of  reasonable  size,  the  DFT  is  implemented  using  an  algorithm 
called  a  fast  Fourier  transform  (FFT).  Older  implementations  of  the  FFT  required  N  to 
be  a  power  of  2;  if  N  wasn’t  a  power  of  2,  zeroes  were  typically  appended  to  the  data  set 
to  increase  its  length  to  a  power  of  2.  Of  course,  this  zero  padding  changed  the  data,  and 
so  had  the  obvious  side  effect  of  introducing  spurious  frequencies  to  the  spectrum.  These 
unwanted  frequencies  were  tolerated  some  years  ago  because  only  by  zero  padding  to  a 
power-of-2  length  could  the  FFT  be  used  at  all.  But  current  FFT  algorithms  are  very 
fast  even  when  N  is  not  a  power  of  2,  and  so  the  absolute  need  for  zero  padding  is  now  a 
thing  of  the  past.  (Appendix  D  has  more  discussion  of  this  zero  padding.)  Note  that  this 
padding  is  distinct  from  the  appending  with  zeroes  that’s  required  when  using  DFTs  to 
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convolve  sequences  of  different  lengths,  discussed  in  Appendix  B,  and  which  has  nothing 
to  do  with  powers  of  2.  That  procedure  works  exactly  with  the  extra  zeroes:  it’s  simply 
not  about  building  a  spectrum  or  a  choice  of  FFT  algorithm. 

Nonetheless,  emitting  a  power-of-2  number  of  pulses  in  a  coherent  processing  interval 
still  seems  to  be  common  practice,  even  in  digital  radars. 

The  Matlab  function  fft  orders  its  DFT  amplitudes  in  a  way  that  renders  negative 
frequencies  positive  and  greater  than  the  Nyquist  frequency.  I  don’t  see  any  point  in 
adding  twice  the  Nyquist  frequency  to  a  negative  frequency  purely  to  make  it  positive  (as 
if  negative  frequencies  are  something  to  avoid),  when  that  resulting  positive  frequency, 
being  now  greater  than  Nyquist,  is  invalid  from  a  DFT  point  of  view.  Instead,  I  give  an 
example  of  implementing  the  DFT  in  Matlab  that  leaves  negative  frequencies  negative, 
using  the  code  of  Table  Cl  in  Appendix  C.  Note  that  the  arrays  x  and  X  in  that  table  are 
arrays  of  numbers;  that  is,  the  first  element  of  X  is  denoted  X(l)  in  Matlab,  even  though 
it’s  called  -X_nv/2i  in  (6-1).  Table  C2  gives  equivalent  Mathematica  code. 

Here  is  an  example  of  plotting  a  frequency  spectrum  using  Table  Cl’s  Dft  function  in 
Matlab.  Suppose  the  data  has  been  sampled  in  time  intervals  of  ts  into  an  array  data. 
The  DFT  of  data  is  the  array 

DFT_of_data  =  Dft(data); 


Plot  DFT_of  _data  versus  frequency  by  first  constructing  an  appropriate  frequency  domain, 
which  will  be  a  sequence  of  N  numbers  {— /,  — /  +  A/,  — /  +  2A /, . . .  }  for  some  /  >  0 
and  frequency  increment  A/.  The  increment  A /  is  the  reciprocal  of  the  duration  of  the 
sampling  interval.  Note  that  this  duration  is  not  the  time  tj  of  the  last  sample  minus  the 
time  tt  of  the  first  sample.  Rather,  the  periodic  nature  of  sinusoids  means  that  the  internals 
of  the  DFT  really  work  on  an  infinite-length  array  made  by  repeating  the  data  endlessly 
before  and  after  the  actual  period  of  sampling.  So  the  DFT  effectively  assumes  the  data 
begins  anew  at  time  t, j  +  ts,  and  thus  the  sampling  interval  has  duration  tf  +  ts  —  t^.  If 
the  data  consists  of  N  samples  spaced  ts  apart,  this  duration  must  equal  Nts,  so 


A/ 


1 


tf  +  ts  —  ti 


1 

Ws- 


(6.3) 


The  value  of  /  must  be  determined  separately  for  each  of  the  cases  of  N  even  or  odd. 


N  even:  For  example,  suppose  N  =  4  so  that  (6.1)  gives  the  transform  as  {X_2,  X_1 ,  X0,Xi}. 
It  follows  that  in  general,  the  frequency  domain  is  {— /, —  A/},  and  this  has 
extent  /  —  A /  —  (— /),  which  we  know  must  equal  ( N  —  1)A/.  Hence 


?  _  NAf  (6.3)  J_ 
2  2ts 


Nyquist  frequency. 


(6.4) 


N  odd:  Suppose  N  =  5,  so  that  the  transform  is  {X_2,  X_1,  X0,  Xi,  X2}.  In  general  the 
frequency  domain  is  {— /, ...,/}  with  extent  2/,  which  must  also  equal  ( N  —  1)A /. 
Hence 


i  (A-l)A/ 
2 


NAf  A  / 
2  2 


Nyquist  frequency 


A / 
2 


(6.5) 


This  can  all  be  implemented  in  Matlab  with 
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Figure  5:  Left:  Data  sampled  at  1^0  Hz  of  two  complex  exponentials  at 
frequencies  5  Hz  and  —12  Hz.  Right:  The  DFT  has  peaks  at  the  correct 
frequencies  with  the  correct  weights. 


Delta_f  =  samplingFrequency/N ; 
nyquistFrequency  =  samplingFrequency/2 ; 

if  abs ( mod (N . 2 ) )  <  le-10 

"/o  DFT  has  even  length  N. 
f_hat  =  nyquistFrequency; 

f requencyDomain  =  linspace(-f  hat,  f _hat -Delta_f ,  N); 
else 

*/.  DFT  has  odd  length  N. 

f_hat  =  nyquistFrequency  -  Delta_f/2; 
f requencyDomain  =  linspace(-f  hat.  f_hat ,  N); 
end 

The  plot  of,  say,  the  absolute  value  of  the  DFT  is  accomplished  with 
stem ( f requencyDomain ,  abs ( DFT_of_data) ) 

As  an  example,  suppose  we  sample  the  function  ei27r(-r>  —  3  e127r(“12  intervals 

of  1/40  second.  The  Nyquist  frequency  is  half  the  sampling  rate,  or  20  Hz,  so  the  DFT 
will  certainly  detect  the  frequencies  of  5  Hz  and  —12  Hz.  Suppose  now,  that  by  chance  we 
happened  to  choose  the  end  sampling  time  to  be  one  sampling  interval  (1/40  s)  less  than 
a  period  of  time  that  contains  a  whole  number  of  cycles  of  both  frequencies  present.  One 
second  is  such  a  period  of  time  (it  contains  exactly  5  cycles  of  one  wave  and  12  cycles  of 
the  other),  so  suppose  we  sample  from  time  0  to  time  1  —  >/40 s,  shown  at  left  in  Figure  5. 
This  sampling  choice  guarantees  that  the  DFT  “sees”  two  infinitely  long  sinusoids,  so  that 
it  will  faithfully  produce  only  peaks  at  the  two  frequencies  5  Hz  and  —12  Hz.  This  can 
indeed  be  seen  in  the  frequency  spectrum  at  the  right  in  Figure  5,  where  we  have  plotted 
the  actual  DFT  and  not  its  absolute  value,  because  in  this  case  it’s  composed  of  real 
numbers  anyway.  The  peaks  are  in  the  expected  places  with  the  correct  weights. 

Now  suppose  we  sample  the  same  function  at  intervals  of  1/20  second,  from  0  to 
1  —  Y20S,  shown  at  left  in  Figure  6.  The  Nyquist  frequency  is  now  10  Hz.  This  is  high 
enough  to  capture  the  5  Hz  sinusoid  (shown  at  the  right  in  Figure  6),  but  not  high  enough 
to  capture  the  one  at  —12  Hz.  As  a  result,  the  DFT  algorithm  effectively  shifts  the 
—  12  Hz  peak  up  or  down  by  multiples  of  the  frequency  domain  width  (which  is  twice  the 
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Figure  6:  Left:  Data  sampled  at  20  Hz  of  two  complex  exponentials  at  fre¬ 
quencies  5  Hz  and  —12  Hz.  Right:  The  lower  sampling  rate  preserves  the  5  Hz 
peak,  but  the  —12  Hz  peak  has  been  aliased  to  —12Hz+  sampling  rate  =  8  Hz. 


Nyquist  frequency — i.e. ,  the  sampling  frequency)  until  the  result  is  within  the  frequency 
domain.  So  we  find  the  correct  frequency  “aliased”  to  —12  Hz  +  20  Hz  =  8  Hz,  as  seen  at 
the  right  in  Figure  6. 

More  realistically,  suppose  now  that  we  sample  the  same  function  from  0  to  fully  1  sec¬ 
ond,  since  we’ll  generally  have  no  knowledge  of  what  frequencies  make  up  the  spectrum— 
that  is,  after  all,  why  we  are  Fourier  transforming  the  data  in  the  first  place!  This  will 
add  one  extra  sampling  point  to  the  sinusoids,  with  the  result  that  the  DFT  no  longer 
sees  two  pure  sinusoids;  rather,  it  sees  a  jump  where  each  copy  of  the  signal  begins  and 
ends,  in  the  infinite  chain  of  copies  placed  side  by  side  in  time  that  the  DFT  works  with. 
As  a  result,  even  if  we  do  sample  quickly  enough  to  capture  both  frequencies — again  at 
1/40  second  at  top  left  in  Figure  7 — the  discontinuity  in  the  joining  of  the  sinusoids  in  the 
data  results  in  many  more  pure  sinusoids  (i.e.  frequencies)  being  present  in  the  Fourier 
spectrum.  That  is,  more  sinusoids  are  required  to  synthesise  a  curve  with  a  discontinuous 
jump. 

This  sampling  behaviour  of  the  discrete  Fourier  transform  is  the  reason  why  more 
frequencies  are  present  (now  with  complex  weights)  in  Figure  7  around  the  dominant 
frequencies  of  5  Hz  and  —12  Hz.  This  broadening  of  the  peaks  in  the  Fourier  spectrum  is 
an  unavoidable  consequence  of  our  practical  inability  to  sample  the  data  set  for  just  the 
right  amount  of  time  to  ensure  each  sinusoid  joins  smoothly  to  a  copy  of  itself  when  the 
data  is  repeated  endlessly  by  the  Fourier  transform. 


6.1  Windowing  Sampled  Data 

The  above  jump  that  occurs  in  the  wrap-around  of  the  sampling  sequence  will  almost 
always  be  present  in  any  real  sample,  because  we  almost  never  deal  with  simple  sinusoids 
along  with  time  intervals  that  have  fortuitously  been  chosen  as  accurately  as  was  done 
to  create  Figures  5  and  6.  The  extraneous  frequencies  introduced  by  this  jump  can  be 
partly  removed  by  eliminating  the  jump;  that  is,  by  forcing  the  data  to  match  up  in 
the  sampling  period  wrap-around  that  the  Fourier  transform  “assumes”  takes  place.  The 
usual  way  of  doing  this  is  simply  to  reduce  the  data  values  to  zero  near  the  beginning 
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Figure  7:  Top  left:  Data  sampled  at  jOHz  of  two  complex  exponentials  at 
frequencies  5  Hz  and  —12  Hz.  Top  right:  The  data  has  been  sampled  in  such 
a  way  as  to  cause  a  jump  where  the  DFT  “joins”  copies  of  it  together  in  an 
infinite  chain.  The  result  is  that  more  frequencies  are  required  to  reconstruct 
the  data.  Bottom:  Absolute  values  of  the  complex  weights  of  the  spectrum 
at  top  right. 


and  end  of  the  sampling  period.  This  is  done  by  multiplying  the  data  by  a  smooth  bell¬ 
shaped  function  that  has  value  zero  at  the  end  points  of  the  data  sequence  and  value 
one  in  its  middle.  Several  choices  of  this  “windowing”  function  exist,  each  with  different 
trade-offs  in  the  resulting  spectrum.  Perhaps  the  two  most  popular  are  the  Hamming 
and  Blackman  windows.  The  three  main  indicators  of  filter  performance  in  the  frequency 
domain  are: 

Flat  passband:  Ideally,  frequencies  required  to  be  passed  with  equal  weightings  should 
not  acquire  any  ripples  in  those  weightings.  The  Blackman  window  has  less  passband 
ripple  than  the  Hamming. 

High  roll-off  rate:  There  should  be  a  well-defined  drop  in  the  Fourier  spectrum  between 
frequencies  required  to  be  passed  and  frequencies  required  to  be  stopped.  The  Ham¬ 
ming  window  has  a  slightly  faster  (better)  roll-off  than  the  Blackman. 

Good  stopband  attenuation:  Frequencies  required  to  be  stopped  by  the  filter  should 
be  attenuated  very  strongly.  The  Blackman  window  has  a  better  stopband  attenua¬ 
tion  than  the  Hamming. 
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Figure  8:  Left:  The  sampled  data  of  Figure  7,  now  Blackman  windowed. 

Right:  The  neighbouring  frequencies  near  5  Hz  and  —12  Hz  have  been  sup¬ 
pressed,  at  the  cost  of  the  peaks  at  5  and  —12  being  widened. 

The  stopband  attenuation  is  arguably  more  important  than  roll-off,  so  the  Blackman 
window  might  be  considered  as  generally  the  most  desirable  to  use,  as  discussed  in  [5]. 

If  we  take  the  “badly  sampled”  data  plotted  in  Figure  7  and  create  a  new  data  set 
by  multiplying  this  data  by  a  Blackman  window,  and  then  plot  the  absolute  value  of  the 
DFT  of  this  new  data,  we  obtain  the  plots  of  Figure  8.  The  windowing  does  suppress  the 
spurious  “side  lobes”  around  5  Hz  and  —12  Hz,  but  at  the  cost  of  broadening  and  lowering 
the  two  peaks  at  those  frequencies. 

Windowing  is  used  in  range-Doppler  processing  to  reduce  side  lobes  in  both  Doppler 
and  range.  For  Doppler  windowing,  recall  that  for  each  range  value,  the  Doppler  infor¬ 
mation  in  the  received  pulse  is  simply  the  Fourier  spectrum  of  the  column  of  correlations 
in  (5.9).  As  discussed  in  the  previous  few  pages  and  in  Figure  8,  we  can  eliminate  spuri¬ 
ous  Doppler  frequencies  by  windowing  each  of  these  columns  before  Fourier-transforming 
them. 

Windowing  in  range  can  be  performed  by  windowing  a  copy  of  the  emitted  pulse,  then 
correlating  that  copy  with  the  returned  pulse.  For  modelling  the  radar  interaction,  this 
procedure  modifies  (4.10)  in  the  following  way.  Write  the  windowed  emitted  pulse  as  Ew, 
in  which  case  (4.10)  becomes 

correlation  =  Ew  *  R  =  E\l.  *  R  =  E\u  *  (E  *  h)  =  {E\,  *  E)  *  h  .  (6.6) 

As  with  the  unwindowed  case,  Ew  *  E  can  be  calculated  offline. 

Range-Doppler  plots  showing  “before  and  after”  results  of  Doppler  windowing  are  given 
in  Section  10.  The  use  of  range  windowing  does  not  give  as  dramatic  an  improvement  in 
the  final  plot  as  the  use  of  Doppler  windowing,  so  no  range  windowing  has  been  used  in 
those  plots. 


7  Representative  Pulse  Types 

A  rectangular  pulse  can  be  viewed  as  a  continuous  wave  being  switched  on  for  a  time  r, 
its  pulse  width.  Aside  from  this,  we  investigate  two  other  pulse  types  in  this  report. 
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7.1  Barker-coded  Pulse 

A  Barker-coded  pulse  has  more  structure  than  a  rectangular  pulse  of  the  same  width. 
Signals  with  rich  structure  autocorrelate  very  well  (i.e.  they  correlate  with  themselves  to 
give  a  well-defined  and  narrow  single  peak,  resulting  in  a  high-quality  matched  filter),  but 
the  larger  bandwidth  required  to  build  this  structure  destroys  some  Doppler  information 
when  compared  to  the  smaller  bandwidth  of  a  rectangular  pulse.  The  reason  is  that  a 
Doppler  frequency  shift  will  generally  be  less  apparent  (i.e.  less  measurable)  when  a  larger 
spread  of  frequencies  is  present. 

The  modulation  that  defines  a  Barker  code  is  a  series  of  rectangular  pulses,  called 
chips,  placed  side  by  side  with  no  gap,  with  weightings  either  +1  or  —1;  that  is,  they 
either  have  no  effect  on  the  phase  of  the  carrier  (+1)  or  they  add  180°  to  it  (—1).  Such  a 

pulse  autocorrelates  to  give  a  high  central  peak  when  the  signal  matches  with  itself,  flanked 

by  low  correlation  “side  lobes”  with  a  saw-tooth  shape.  Demanding  an  even  saw-tooth 
pattern  of  side  lobes  results  in  just  seven  known  Barker  codes:  it  has  been  conjectured 
that  no  longer  Barker  codes  exist.  For  an  example  of  Barker  autocorrelation,  consider  the 
length-7  Barker  code:  [+1  +1  +1  —1  —1  +1  —  l]  (i.e.  this  has  7  chips).  This  autocorrelates 
to  give  the  following,  remembering  from  Appendix  B  that  correlation  is  accomplished  by 
complex-conjugating  one  reversed  signal  (which  is  real  in  this  case)  and  convolving  this 
with  the  other  signal — which  here  is  of  course  identical  to  the  first  signal: 

[-1  1-1-1  1  1  1]  *  [1  1  1-1-1  1  -1] 

=  [-1  0  -1  0  -1  0  7  0  -1  0  -1  0  -1]  .  (7.1) 

We  see  here  the  high  peak  in  the  middle,  with  the  characteristic  Barker  saw-tooth  side 
lobes. 


7.2  Chirped  Pulse 


A  chirped  signal  varies  its  frequency  continuously  throughout  each  of  its  pulses.  This 
endows  it  with  much  structure,  resulting  in  a  higher-quality  autocorrelation.  The  sin¬ 
gle  chirped  pulse  thus  acquires  the  range  resolution  that  only  a  much  shorter  constant- 
frequency  pulse  would  have.  But  the  benefit  of  using  the  longer  (chirped)  pulse  is  that 
this  long  pulse  requires  a  lower  energy  density  than  a  short  (non-chirped)  pulse  to  achieve 
the  same  result,  and  so  is  both  easier  to  generate  and  less  visible  to  the  target  than  the 
non-chirped  pulse. 

As  in  (2.13),  write  the  emitted  signal  as  sc(t)  =  u(t )  elu}°t.  We  will  show  shortly  that 
for  a  signal  beginning  at  time  t0  composed  of  N  chirped  pulses  of  period  T  and  width  r, 
the  emitted  chirped  signal  is 


sc(t)  =  e 


i  LJ0t 


F(t  —  to)  +  F(t  —  to  —  T)  +  •  •  •  +  F(t  —  to  —  N  —  1 T)  , 


(7.2) 


where 


m 


gi^t2  O^t^T 

0  otherwise 


(7.3) 


for  some  constant  /r.  (In  fact,  it  will  turn  out  that  such  an  sc(t)  is  not  quite  an  analytic 
signal  because  its  spectrum  can  contain  negative  frequencies;  however,  these  are  negligibly 
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present  when  the  carrier  frequency  is  large  enough,  which  it  usually  is.  See  further  discus¬ 
sion  and  examples  in  Section  9.1.)  The  rationale  behind  the  form  of  the  modulation  F(t) 
in  (7.3)  is  as  follows.  A  general  signal  sc(t)  =  Aft )  can  be  Fourier  decomposed  as 

OO 

scif)  =  J  Sc(u)  exut  do; ,  (7.4) 

—  OO 

where 

OO  OO 

Sc(u)  =  ^  J  sc(t )  e-^  dt  =  ±  j  A(t)  dt  (7.5) 

— OO  — OO 

Now  choose  any  value  of  u:  almost  no  contribution  to  the  last  integral  in  (7.5)  occurs 
when  its  phase  (j>(t)  —  ut  is  changing  rapidly  with  time,  because  in  such  a  case  many 
complex  numbers  of  essentially  random  phases  are  being  added,  giving  a  result  of  about 
zero.  This  cancellation  won’t  happen  when  <f(t)  —  ut  is  not  changing  rapidly  with  time. 
That  is,  we  expect  the  dominant  contribution  to  the  integral  to  occur  when  (f>(t )  —  ut  is 
roughly  constant  over  time,  or  d/d t  [(/(f)  —  ut]  ~  0,  or  (ft if)  ~  u.  (This  reasoning  is  known 
as  the  method  of  stationary  phase.)  For  example,  at  some  time,  say  t  =  41,  Sc(u)  has  the 
most  support  at  u  =  (/'(41)  =  (say)  2.  That  is,  we  presume  that  at  t  =  41,  the  frequency 
spectrum  Sc(u)  is  peaked  around  a  dominant  frequency  of  u  =  2.  In  general,  at  time  t 
the  dominant  frequency  of  the  signal  sc(t)  =  is 


(7.6) 


Refer  to  (7.2)  and  (7.3)  and  focus  on  the  basic  signal  component  At  time  t  its 

dominant  angular  frequency  is  u0  +  2nfj.t,  which  is  the  carrier  plus  an  extra  part  that  grows 
linearly  with  time;  this  constitutes  the  chirp  and  is  the  rationale  for  the  form  of  (7.3).  The 
angular  frequency  u0  +  2irfit  corresponds  to  an  actual  frequency  of  /0  +  fat,  which  explains 
why  7 r  is  conventionally  included  in  (7.3).  It’s  clear  that  fi  is  the  slope  of  the  straight-line 
plot  of  chirped  frequency  versus  time,  as  shown  in  Figure  10  on  page  36. 


8  Selecting  Valid  Pulse  Parameters 

In  this  section  we  describe  how  to  choose  pulse  parameters  that  guarantee  valid  range- 
Doppler  plots.  We  illuminate  a  single  possibly  moving  target  at  distance  r  from  our  radar 
by  the  train  of  pulses  shown  in  Figure  9.  The  parameters  describing  this  pulse  train  are 
used  in  the  subsections  that  follow. 

8.1  Attainable  Range  Measurement  of  a  Pulse  Train 

A  radar’s  receiver  is  necessarily  switched  off  during  the  time  interval  r  that  it  emits  a  pulse. 
Any  target  closer  than  cr/2  will  then  cause  some  of  the  reflected  pulse  to  be  received  while 
the  receiver  is  switched  off.  This  minimum  range  cr/2  is  called  the  radar’s  blind  range. 

Likewise,  the  return  from  any  target  at  range  greater  than  cT/2  (where  T  is  the  PRI) 
will  arrive  at  the  radar  receiver  either  during  or  after  it  emits  the  next  pulse  of  its  train. 
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Pulse  width 


Pulse  rep.  interval  T 


sampling  interval  ts 


time 


Figure  9:  A  series  of  pulses  with  parameters  indicated.  Only  the  envelopes 
of  the  pulses  are  shown. 

We  cannot  know  that  such  a  return  belonged  to  the  first  pulse  and  not  the  second,  so 
ranges  greater  than  cT / 2  create  ambiguities  in  the  measured  range;  thus  cT / 2  is  called 
the  maximum  unambiguous  range.  The  radar’s  useful  range  extent  then  runs  from  the 
blind  range  rmin  =  cr/2  to  the  maximum  unambiguous  range  rmax  =  cT / 2: 


(8.1) 


Note  that  rmin  =  cr/2  is  also  the  range  resolution  or  range  bin  width  A rres  of  a  simple 
rectangular  pulse  [i.e.  p  =  0  in  (7.3)],  being  the  difference  of  the  two  closest-together 
ranges  r\ ,  r2  such  that  the  returned  pulse  from  a  target  at  r1  will  be  received  without 
overlapping  the  returned  pulse  from  a  target  at  r2.  Section  9 — starting  with  (9.1) — 
calculates  the  range  resolution  of  pulses  with  more  structure. 


cr/2  ^  r  ^  cT /2  . 


8.2  Bounds  on  the  Pulse  Repetition  Interval  T 

Given  a  target  at  range  r,  the  above  discussion  of  blind  range  and  maximum  unambiguous 
range  shows  that  r  and  T  must  be  chosen  to  ensure  the  target  is  visible  to  the  radar,  or 

ct/2  ^r  ^cT/2.  (8.2) 

This  places  bounds  on  the  pulse  width  r  and  the  PRI  T : 

crmax/2  =  r  =  cTmin/2 .  (8.3) 

We  focus  in  particular  on  Tmin: 

Tmin  =  2r/c.  (8.4) 

An  upper  bound  Tmax  to  the  PRI  arises  by  avoiding  Doppler  ambiguities  in  the  following 
way.  Refer  to  (2.27),  which  gives  the  Doppler  frequency  as 

/d  = -2/ ov/c  (8.5) 
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for  a  carrier  frequency  /0.  We  will  be  extracting  the  Doppler  frequency  from  numbers 
sampled  at  a  rate  of  1/T  ( not  1  /ts:  remember  that  Doppler  information  is  sampled  from 
pulse  to  pulse,  at  one  range  value).  Nyquist’s  theorem  says  that  approximately  half  this 
rate  is  the  maximum  Doppler  frequency  that  a  DFT  can  detect,8  so 

(8.6) 


This  rearranges  to  give 


T  < 


c 

4/b  I  ^ 


so  that  Tmax 


c 

4/oh 


We  must  choose  a  PRI  between  Tmin  and  Tmax: 


(8.7) 


2  r 
c 


<  T  < 


c 

4/oh 


(8.8) 


8.3  Attainable  Speed  Measurement  of  a  Pulse  Train 


There  is  a  maximum  unambiguous  speed  umax  that  a  given  pulse  train  is  able  to  detect, 
irrespective  of  whether  the  target  is  approaching  or  receding.  Also,  while  the  train  can  in 
principle  detect  a  minimum  speed  of  zero,  it  does  have  a  speed  resolution  Aures.  We  will 
calculate  umax  and  Aures  here.  First,  (8.5)  gives  the  target’s  radial  velocity  (positive  for 
receding)  as 


v  = 


-c/p 
2/o  ' 


(8.9) 


Remember  from  (8.6)  that  the  maximum-detectable  Doppler  frequency  is  approximately 
/ma,x  =  l/(2T).  The  maximum  unambiguous  speed  comes  from  the  absolute  value  of  (8.9): 


sr 


c/ 


max 

D 


2/0 


cl  c 
2/o  2T  =  4/gT  ' 


(8.10) 


The  speed  resolution  Aures  arises  because  the  calculated  Doppler  spectrum  is  a  set  of 
points  and  not  a  continuous  curve,  since  this  spectrum  results  from  Fourier-transforming 
the  N  numbers  in  (5.9).  The  resolution  corresponds  to  the  speed  spacing  between  these 
points.  The  Doppler  frequency  spacing  is  A  fD  =  1  /(NT)  [use  (6.3)  with  ts  replaced  by  T ], 
and  this  converts  to  the  speed  resolution  using  the  same  factor  as  in  (8.9): 


Aures  = 

-cA/d 

c  c 

2/o 

2 foNT  2/0AfCPI  ’ 

(8.11) 


where  AfCPI  is  the  length  NT  of  the  CPI. 


8I  say  “approximately”  because  the  actual  number  is  /  in  (6.4)  and  (6.5),  where  the  ts  of  those  equations 
is  replaced  by  T  here. 
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Alternative  Views  of  umax  and  Avres 


It’s  instructive  to  derive  the  maximum  unambiguous  speed  detectable  by  one  particularly 
simple  wave  train  in  an  alternative  way  to  the  above.  Suppose  the  radar  emits  an  infinite 
number  of  impulses  spaced  T  apart.  (That  is,  let  t  — >•  0  and  N  — >•  oo.)  The  frequency 
spectrum  of  this  train  is  easily  calculated:  it  turns  out  to  be  an  infinite  series  of  spikes 
spaced  1/T  apart.  We  can  define  the  maximum  unambiguous  speed  to  be  that  whose  effect 
is  to  shift  these  spikes  by  no  more  than  half  their  spacing,  since  otherwise  we  wouldn’t  be 
able  to  tell  whether  any  particular  spike  had  moved  to  the  left  or  the  right.  That  is,  we 
require  the  absolute  value  of  the  Doppler  shift,  f02v/c  from  (2.27),  to  be  less  than  half 
the  frequency  spike  spacing  of  1/T: 


fp2v  ^  1 
c  "  2T  ' 


(8.12) 


Equality  in  (8.12)  defines  vmax: 


max  _ 

~~ c  =  2  T  ’ 


(8.13) 


and  this  rearranges  to  give  (8.10)  again. 

The  speed  resolution  Aures  is  sometimes  related  heuristically  to  a  simple  wave  picture  in 
the  following  way.  (This  argument  doesn’t  use  the  above  infinite  wave  train.)  Recall  (8.11) : 


AUes 


cA/d 

2/o 


(8.14) 


where  A  fjj  is  the  (positive)  spacing  between  points  in  the  Doppler  spectrum — which  is 
also  the  smallest  positive  Doppler  frequency  measurable.  Let  this  smallest  frequency 
correspond  to  a  nominal  longest  “Doppler  wavelength”  XD  =  c/AfD  (but  note  that  this 
does  not  correspond  to  a  physical  wave).  Because  A fD  =  1  /{NT)  =  l/AfCPI,  we  have 


\D  — 


c 


c  AfCPI . 


(8.15) 


This  equation  portrays  the  speed  resolution  as  resulting  from  the  fitting  of  precisely  one 
of  these  longest  “Doppler  wavelengths”  into  the  distance  travelled  by  light  in  one  coherent 
processing  interval. 


DistanceAAlocity  Maxima  Tradeoff 

Equations  (8.1)  and  (8.10)  combine  to  give 


T  V 
'  max  17  max 


(8.16) 


We  see  here  a  tradeoff  between  maximum  unambiguous  range  and  maximum  unambigu¬ 
ous  speed:  for  a  given  carrier  frequency  /0,  increasing  one  must  decrease  the  other.  This 
is  another  reason  why  high  carrier  frequencies  are  desirable  in  radar:  a  high  carrier  fre¬ 
quency  /o  allows  both  rmax  and  umax  to  be  large  simultaneously.  But  see  the  remark  just 
after  (8.18)  ahead. 
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8.4  Doppler  Aliasing 

The  target’s  radial  velocity  was  given  in  terms  of  the  Doppler  frequency  in  (8.9).  But  all 
velocities  that  produce  Doppler  frequencies  outside  the  DFT’s  domain  will  be  aliased  into 
that  domain  by  the  action  of  the  DFT:  the  DFT  will  add  the  necessary  number  of  multiples 
of  the  frequency  domain  width  1/T  to  the  Doppler  frequencies  of  those  higher  velocities 
to  shift  their  frequencies  into  that  domain.  In  other  words,  the  Doppler  frequency  f°FT 
produced  by  the  DFT  (which  might  not  equal  the  actual  Doppler  frequency  fjj)  will  be 
indistinguishable  from  all  frequencies  present  in  the  data  that  have  value  f°FT  +  n/T  with 
n  ranging  over  all  integers.  Those  other  frequencies  correspond  to  velocities  vn  via  (8.9): 


v 


n 


-c(/DDFT  +  n/T) 

2/o 


for  all  integers  n. 


(8.17) 


So  all  vn  will  be  aliased  to  the  v  given  by  (8.9),  meaning  that  a  target  with  recession 
velocity  vn  (for  any  integer  n)  will  appear  to  have  velocity  v  on  the  range-Doppler  plot.  In 
particular,  target  velocities  that  are  aliased  down  to  zero  speed  are  potentially  hazardous, 
since  they  become  indistinguishable  from  the  background  and  hence  cannot  trigger  any 
velocity-sensitive  sensor  to  acquire  the  target  as  being  possibly  of  interest.  These  blind 
speeds  are  |un|  such  that  f°FT  =  0,  or 


CTl 

\vn\  ^2/oT’  for  n=  (8-18) 

Unfortunately,  given  that  the  spacing  between  blind  speeds  is  c/(2/0T),  a  high  carrier 
frequency  /0  allows  more  blind  speeds  in  some  given  speed  range.  This  is  a  counter¬ 
argument  to  the  remark  just  after  (8.16)  that  high  carrier  frequencies  are  desirable. 


8.5  Range— Doppler  Coupling 

Refer  to  Figure  10,  whose  black  line  segment  is  the  frequency  versus  time  of  a  single  chirped 
pulse  of  duration  r,  the  pulse  width.  A  stationary  target  returns  a  pulse  whose  frequency 
is  plotted  as  the  blue  line  segment,  and  a  moving  target  at  the  same  range  returns  a  pulse 
whose  frequency  is  plotted  as  the  red  line  segment.  The  amount  of  Doppler  shift  fD  in  the 
picture  is  exaggerated;  a  realistic  figure  is  about  one  millionth  of  the  carrier  frequency  /0. 

The  slope  of  the  red  segment  is  actually  slightly  different  to  that  of  the  blue:  if  the  frequency  of 
the  emitted  pulse  is  /  =  f0  +  /rt  (refer  Section  7.2),  and  the  returned  pulses  start  at  time  td  (the 
time  delay  in  Figure  1),  then  the  blue  segment  is  just  the  black  segment  shifted  to  the  right  by  td, 
meaning  the  blue  has  equation 

f  =  fo+ ■  (8.19) 

This  has  slope  fi-  But  (2.25)  says  that  the  received  frequency  is  o>0(l  —  2v/c),  so  that  the  red 
segment  has  equation 

f  =  (1  -  2 v/c)  [fo  +  -  td)] ,  (8.20) 

which  has  slope  (1  —  2v/c)n-  The  value  of  v/c  is  typically  ±10~6,  so  the  slopes  of  red  and  blue 
are  almost  identical. 

The  key  feature  of  Figure  10  relates  to  finding  the  parts  of  the  emitted  and  received 
pulses  that  correlate  best.  The  emitted  pulse  segment  labelled  “1”  in  the  figure  (the  small 
black  rectangle)  will  correlate  well  with  the  returned  segment  labelled  2  in  the  absence  of 
Doppler  shift.  In  the  presence  of  Doppler  shift,  we  might  hope  segment  1  would  correlate 
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Figure  10:  Range-Doppler  coupling  for  a  chirped  pulse.  The  black  line 
segment  is  the  emitted  pulse.  A  returned  pulse  with  no  Doppler  shift  is  the 
blue  line  segment,  and  a  returned  pulse  with  Doppler  shift  is  the  red  line 
segment,  in  this  case  for  an  approaching  target. 


well  with  segment  2l  in  order  to  derive  the  same  time  shift  necessary  to  give  the  range. 
But,  in  fact,  segment  1  might  better  correlate  with  segment  3  since  they  share  the  same 
frequency.  (This  assumes  some  chirp  is  present,  in  which  case  we  cannot  expect  to  be  able 
to  take  the  limit  g  — >  0  at  the  end  of  our  calculation.)  That  means  the  Doppler  shift  has 
the  effect  of  making  it  appear  as  if  the  returned  pulse  arrived  early  by  an  error  time  t E,  for 
the  case  when  fp  has  the  same  sign  as  g.  (If  the  signs  are  different,  the  red  segment  will 
appear  to  have  arrived  later  than  the  blue  segment.)  Because  the  red  line  in  the  figure  has 
about  the  same  slope  as  the  blue  and  black  lines,  we  can  write  g  ~  /dAb  for  any  choice 
of  signs,  and  infer  a  range  error  of 


A r  =  measured  range  —  true  range 
=  c(td  -  tE)  _  eta 
2  2 

=  ~ctE  =  ~cfD  =zflx  ~2vfo  =  YA 
2  2g  2  g  c  g 


(8.21) 


This  error  is  called  range-Doppler  coupling ,  and  amounts  to  about  7  centimetres  for  the 
target  of  Figure  16  ahead.  Given  that  it  has  such  a  small  value  for  our  scenarios  of  interest, 
we  don’t  consider  it  further  in  this  report. 


9  Range  Resolutions  of  Different  Pulse  Types 

A  standard  expression  in  radar  signal  processing  is  that  a  pulse’s  range  resolution,  or  range 
bin  width,  is  given  by 

Arres  ~  ^  ,  (9.1) 

where  B  is  the  pulse’s  bandwidth.  This  isn’t  a  hard  and  fast  rule,  in  that  it  suggests 
that  “bandwidth”  is  a  well-defined  term,  which  isn’t  quite  the  case.  In  principle  the  total 
bandwidth  of  a  signal  is  defined  as  the  width  (say,  the  “full  width  at  half  maximum”)  of  the 
support  of  the  signal’s  frequency  spectrum;  but  this  width  can  be  infinite,  so  in  practice  the 
total  bandwidth  is  defined  as  the  width  of  the  dominant  part  of  the  spectrum’s  support. 
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One  rationale  for  (9.1)  is  the  following.  Most  pulses  used  in  radar  have  a  spectrum 
that  is  shaped  something  like  a  top  hat  of  width  B,  meaning  this  spectrum’s  support  is 
fairly  well  defined  and  the  spectrum  is  reasonably  flat  throughout  this  width.  Neglecting 
noise,  the  returned  pulse’s  spectrum  is  similar  to  the  emitted  pulse’s  spectrum,  save  for  a 
small  Doppler  shift.  The  correlation  of  the  emitted  signal  E  with  the  returned  signal  R 
is  E  -k  R  =  E^  *  R;  in  other  words,  a  convolution  of  a  reversed  and  conjugated  E  with  R. 
This  convolution  is  the  inverse  Fourier  transform  of  a  product:  the  product  of  the  Fourier 
transform  of  and  the  Fourier  transform  of  R,  probably  with  some  zero  appending 
as  discussed  around  (Bll),  which  will  alter  the  bandwidths  somewhat,  but  presumably 
not  by  much.  The  Fourier  transform  of  E'  has  bandwidth  B,  and  the  Fourier  transform 
of  R  has  bandwidth  approximately  B,  and  both  spectra  are  centred  at  almost  the  same 
frequency.  So  their  product  is  approximately  another  top  hat  of  width  B.  The  inverse 
Fourier  transform  of  this  bandwidth- 1?  function  (which  will  be  the  correlation  E-kR)  must 
then  have  a  typical  width  of  l/B.  But  a  well-crafted  pulse’s  correlation  with  its  return 
will  be  a  spike,  so  that  this  spike  has  a  width  of  about  l/B ,  meaning  the  correlation  peak 
is  spread  over  a  time  interval  of  l/B.  The  associated  range  resolution  of  the  pulse  is  c/2 
times  this  width,  which  gives  (9.1). 

Equation  (9.1)  is  reasonable  qualitatively,  in  that  a  large  bandwidth  (large  B )  means 
high  pulse  complexity,  which  means  high  autocorrelation  ability,  which  leads  to  very  good 
range  resolution  (small  Arres).  Fourier  analysis  shows  that  bandwidth  of  a  rectangular 
pulse  formed  by  multiplying  a  single-frequency  carrier  by  a  top-hat  window  of  duration  r  is 
B  ss  1/t,  for  which  case  (9.1)  produces  the  expected  Arres  =  cr/2  quoted  just  after  (8.1). 

So  range  resolution  is  determined  by  bandwidth,  not  pulse  width.  Suppose  we  have 
a  radar  pulse  of  duration  r  and  bandwidth  B,  which  might  have  some  internal  structure 
such  as  a  chirp  or  a  Barker  code.  The  bandwidth  of  this  highly  structured  pulse  is  greater 
than  the  bandwidth  of  a  duration-r  rectangular  pulse  (i.e.  one  with  no  internal  structure 
apart  from  its  carrier),  this  bandwidth  being  1/t.  Now  imagine  a  rectangular  pulse  with 
no  structure,  duration  l/B ,  and  bandwidth  B,  which  has  therefore  the  same  range  reso¬ 
lution  as  the  actual  pulse.  If  the  actual  pulse  has  high  structure  (large  B ),  this  imagined 
rectangular  pulse  with  the  same  large  bandwidth  B  will  have  a  very  small  duration  l/B, 
smaller  than  r.  This  smaller  duration  of  this  (structureless)  rectangular  pulse  is  called 
the  effective  pulse  width  of  the  original  highly  structured  pulse:  reff  =  l/B  <  r.  Although 
this  imagined  narrow  and  simple-to-generate  pulse  has  the  same  range  resolution  as  the 
wide  and  difficult-to-generate  structured  pulse,  emitting  the  imagined  simple  pulse  over 
the  comparatively  short  time  reg  requires  the  radar  to  reach  a  high  power  level,  which 
might  be  impossible;  better  is  for  the  radar  to  emit  the  structured  pulse  over  the  longer 
time  r,  for  which  a  lower  average  power  suffices.  Also,  a  lower  average  emitted  power 
means  the  radar  is  less  visible  to  the  target. 

The  actual  pulse  of  duration  r  is  sometimes  called  the  uncompressed  pulse,  while  the 
imagined  plain  rectangular  pulse  of  duration  reg  is  sometimes  called  the  compressed  pulse. 
The  compression  ratio  is  defined  as  the  ratio  of  these  widths: 

T 

compression  ratio  =  —  =  Bt  ^  1 .  (9-2) 

Aft 


The  compression  ratio  is  a  measure  of  the  emitted  pulse’s  structure:  higher  structure  means 
higher  compression  ratio.  Also,  the  compression  ratio  determines  the  emitted  pulse’s  range 
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resolution  using  (9.1): 


Ar,„  »  4  =  - .  (9,3) 

2 B  2  compression  ratio 

In  words,  the  range  resolution  Arres  of  the  actual  long  structured  pulse  (the  “uncompressed 
pulse”,  duration  r)  equals  the  range  resolution  of  the  imagined  short  structureless  rect¬ 
angular  “compressed  pulse”  (duration  Teff  <  r),  which  equals  the  range  resolution  of  an 
imagined  structureless  rectangular  pulse  (duration  r,  not  reff)  divided  by  the  compression 
ratio. 

In  some  analogue  radars,  chirped  pulses  are  physically  compressed  on  reception  by 
filtering  them  through,  for  example,  a  “surface  acoustic  wave  device”,  which  passes  dif¬ 
ferent  frequencies  at  different  speeds  to  physically  compress  the  returned  pulse,  achieving 
the  range  resolution  of  this  “compressed  pulse”  without  the  high  power  requirement  of 
a  short  pulse.  But  modern  digital  radars  no  longer  compress  pulses  on  reception,  so  for 
them  the  terms  “uncompressed  pulse”,  “compressed  pulse”,  and  “compression  ratio”  are 
not  literal:  the  terms  are  not  meant  to  imply  that  the  received  (or  emitted)  radar  pulse 
is  ever  compressed.  (I  think  this  terminology  of  compression,  once  useful  for  analogue 
radars,  is  misleading  for  digital  radars  and  shouldn’t  be  used.)  Note  also  that  in  some 
texts,  the  actual  pulse’s  autocorrelation  function  is  also  called  the  compressed  pulse,  as 
might  be  the  actual  pulse’s  correlation  with  its  return.  So,  what  is  called  a  compressed 
pulse  might  be  (a)  an  imagined  narrow  pulse  of  duration  rcg,  or  (b)  a  correlation  sequence 
of  the  actual  (“uncompressed”)  pulse  with  either  itself  or  its  return;  and  this  sequence  is 
not  a  radar  pulse  at  all. 

Either  way,  in  a  digital  radar,  a  compressed  pulse  is  not  the  compressed  version  of 
some  kind  of  “uncompressed”  pulse.  Nothing  ever  gets  compressed.  In  the  digital  radar 
age  where  surface  acoustic  wave  devices  are  no  longer  used  to  physically  compress  radar 
pulses,  the  technique  of  frequency  modulation  is  not  meant  to  be  visualised  as  having  the 
magical  effect  of  physically  compressing  a  radar  pulse.  It  simply  gives  the  pulse  more 
structure,  which  leads  to  better  autocorrelation.  The  emitted  or  “uncompressed”  pulse  is 
the  only  radar  pulse  that  actually  exists. 


9.1  Range  Resolution  of  a  Chirped  Signal 

Using  a  large  chirp  [a  large  p  in  (7.3)]  endows  a  pulse  with  more  structure,  which  in  turn 
gives  it  more  bandwidth  (a  “high  compression  ratio”),  thus  giving  a  better  (i.e.  smaller) 
range  resolution. 

Let’s  calculate  the  amount  of  chirp  p  in  (7.3)  required  for  a  given  pulse  width  and 
desired  range  resolution.  The  range  resolution  is  given  by  (9.1): 


A  r 
— res 


C 

2^’ 


(9.4) 


where  we’ve  now  written  Btot  for  the  total  bandwidth  of  a  single  chirped  pulse;  this  band¬ 
width  is  distinct  from  the  amount  by  which  the  pulse’s  instantaneous  frequency  changes 
during  the  chirp,  which  is  its  swept  bandwidth  Bswept.  Although  the  notion  of  bandwidth 
is  not  scrupulously  well-defined,  we  can  find  an  approximate  value  of  this  total  bandwidth 
by  Fourier  transforming  a  single  chirped  pulse,  described  next. 
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Calculating  the  Bandwidth  of  a  Chirped  Pulse 

Over  the  duration  r  of  a  single  pulse,  the  chirped  wave  train  of  (7.2)  and  (7.3)  is 

sc(t)  =  eMfot+^t2)  =  , 


where  /0  is  the  carrier  frequency  and  fj,  is  the  chirp  parameter,  quantifying  how  the  sig¬ 
nal  frequency  changes  during  the  pulse.  Equation  (7.6)  expresses  the  frequency  used  to 
quantify  the  bandwidth  as 

/(*)  =  =  fo  +  l*t  >  (9-6) 

showing  that  the  instantaneous  frequency  f(t)  does  indeed  grow  linearly.  The  pulse’s 
frequency  spectrum  Sc(u)  is  given  by  the  Fourier  transform  of  sc(t): 

/OO  PT  „ 

e-^s.it)  d t  =  /  e^+^)t+^t  d t ,  (9.7) 

-oo  J  0 

where  u>$  =  2irf0.  Evaluate  this  integral  with  the  identity 

/e"“’+to  =  (9-8) 

to  give 


Sc(u) 


erf 


t  + 


i(^  -  o;0) 

2 \J—  Vt Tjl 


erf 


l(u;  -  o;0) 
2  yj—  \1 TfJ, ' 


(9.9) 


Note  that  n  is  not  trivially  able  to  be  set  to  zero  here,  because  the  identity  (9.8)  relies 
on  the  process  of  “completing  the  square”,  which  has  no  meaning  if  /i  =  0.  But  we  can 
certainly  set  /r  to  be  arbitrarily  small. 

Alternatively,  we  could  make  the  real  and  imaginary  parts  of  Sc{u)  more  explicit  by 
using  the  Fresnel  integrals  S(x)  and  C(x),  where 


S(x) 

C(x) 


rx  (  sin  I  7ru2 


o  l cos  J  2 


d  u . 


These  lead  to  the  identity 


where 


eiax2+ibx  dx  =  JL  e^~  [C(U)  +  iS(u)] 

V  ZjCL 


/2a 1  b 
u  =  \  l  —  x  + 


vr  y/2ira'  ’ 

enabling  the  Fourier  transform  in  (9.7)  to  be  written  alternatively  to  (9.9)  as 


5c(w)  = 


]_  —  i(cj— u>p) 


0  47T/X 


+ 


(9.10) 

(9.11) 

(9.12) 


(9.13) 


Plots  of  |5'c(w)|  versus  frequency  /  =  uj/(2tt)  using  (9.9)  are  shown  in  Figure  11  on  the 
next  page.  The  carrier  co0  just  shifts  the  spectrum,  so  has  been  set  to  zero  in  each  plot. 
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Frequency  f 

H  =  2,  r  =  5.  Eqn  (9.17)  gives  Btot  ~  10.2. 


p  =  20,  r  =  5.  Eqn  (9.17)  gives  Btot  ~  100.2. 


Figure  11:  Frequency  spectra  of  single  pulses  with  various  values  of  p  and  t, 
whose  units  are  immaterial  here.  The  carrier  has  been  removed;  thus  zero 
frequency  corresponds  to  the  carrier  frequency.  The  values  of  the  bandwidth 
calculated  using  (9.17)  are  listed,  and  compare  well  with  a  visual  inspection 
of,  say,  the  full  width  at  half  maximum  in  each  plot. 


Note  that  usually  the  bandwidth  is  much  smaller  than  the  carrier  frequency,  so  that  in  accordance 
with  the  comments  in  Section  2  and  Appendix  A,  the  signal  sc(t)  in  (9.5)  essentially  has  only  a 
positive-frequency  spectrum.  Nevertheless,  were  uj0  to  be  very  small,  sc(t)  would  certainly  have 
some  negative  frequencies  in  its  spectrum,  which  is  at  odds  with  the  fact  that  an  analytic  signal 
only  has  positive-frequency  components.  The  fact  is  that  sc(t)  in  (9.5)  is  not  exactly  an  analytic 
signal;  applying  the  Hilbert  transform  to  its  real  part  is  not  guaranteed  to  give  its  imaginary  part. 
But  when  the  carrier  frequency  is  high  enough — which  is  perhaps  always  the  case — the  plots  of 
Figure  11  show  that  the  spectrum  essentially  does  have  only  positive  frequencies. 
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/  =  /o  /  =  /o  +  Vt  f  =  fo  +  HT 


0  T 


Figure  12:  An  idealised  fusion  of  two  plane  wave  trains  used  to  derive  (9.17) 


The  plots  of  Figure  11  show  that  the  behaviour  of  the  spectra  as  g  and  r  vary  is  not 
straightforward.  But  there  is  some  latitude  in  the  definition  of  a  spectrum’s  bandwidth, 
given  that  technically  each  plot  in  Figure  11  has  infinite  support.  The  bandwidth  in  each 
plot  can  be  considered  or  defined  as  a  kind  of  “full  width  at  half  maximum”.  Here  we  derive 
a  simple  expression  for  the  bandwidth  which  agrees  well  with  full-width-at-half-maximum 
values  in  Figure  11. 

Begin  by  treating  a  chirped  pulse  as  the  product  of  a  particular  type  of  idealised  wave 
train  and  a  top-hat  function.  The  idealised  wave  train  is  shown  in  Figure  12.  It  has  a 
constant  amplitude  and  the  following  frequency  behaviour: 

-  For  t  =  —  oo  — >•  0  it  has  a  frequency  f0. 

-  For  f  =  0  — '?  t  it  has  a  frequency  /0  +  gt,  so  that  its  frequency  increases  linearly  at 
a  rate  of  g  over  the  width  r  of  the  pulse,  reaching  a  value  of  f0  +  gr. 

-  For  t  =  t  oo  it  remains  at  this  frequency  of  /0  +  gr. 

The  top-hat  function  is  nonzero  only  for  the  pulse  duration  t  =  0  — >•  r.  The  real  chirped 
pulse  is  then  the  product: 

chirped  pulse  =  idealised  wave  train  x  top-hat  of  width  r.  (9-14) 

The  bandwidth  Btot  is  found  by  Fourier  analysing  this  pulse;  here  ff  {}  stands  for  a  Fourier 
transform  and  for  convolution: 

spectrum  =  ^{chirped  pulse} 

bandwidth  Btot 

=  ^{idealised  wave  train  x  top  hat} 

=  ^{idealised  wave  train}  *  fF { top  hat}  .  (9.15) 

i _ i  i _ i 

bandwidth  =Bsv/ept  a  sine  function 

of  bandwidth  1/r 

The  idealised  wave  train  is  a  fusion  of  two  plane  waves:  its  early-time  or  “left  part” 
has  a  single  frequency  /0  (shown  in  red  in  Figure  12),  its  late-time  or  “right  part”  has 
a  single  frequency  fo  +  gr  (shown  in  blue  in  Figure  12),  and  there  is  an  intermediate 
transition  region  (green)  which  passes  through  all  frequencies  in  between,  f0  — >•  /0  +  gr. 
Remember  that  an  infinite  plane  wave  with  a  single  frequency  has  no  bandwidth,  so  we 
can  approximate  (even  define)  the  bandwidth  of  the  entire  wave  train  as  being  solely  the 
sweep  of  frequencies  in  the  intermediate  green  region: 

£swept  =  Mr-  (9-16) 
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This  is  the  first  factor  in  the  convolution  on  the  right-hand  side  of  (9.15).  The  second  factor 
there  is  a  sine  function  of  well-known  bandwidth  1/r.  Now,  if  we  convolve  two  functions 
of  approximate  widths  \fi\ r  and  1/r  in  (9.15),  we  can  expect  the  result — which  amounts 
to  “smudging”  one  using  the  other9 — to  have  a  width  that  is  approximately  the  sum  of 
these  widths,  or  |/x|r  +  1/t.  So  the  required  single-pulse  bandwidth  Btot  is  approximately 
this  sum: 


Btot  -  Hr  +  1/r  . 


(9.17) 


As  far  as  I’m  aware,  radar  textbooks  tend  not  to  distinguish  between  Bswept  and  Btot, 
and  seem  only  ever  to  write  (9.16)  as  the  expression  for  “the”  bandwidth  B.  But  simply 
writing  “ B  =  |/x| r”  (or  what  is  more  likely,  UB  =  /ur”)  doesn’t  return  the  baseline  value 
of  1/r  in  the  no-chirp  case  (/x  =  0).  In  contrast,  (9.17)  does  return  1/r  in  that  case. 

Applying  (9.17)  to  the  /x,  r  values  in  Figure  11  gives  the  bandwidths  listed  in  its  sub¬ 
captions.  There  is  a  good  visual  agreement  between  (9.17)’s  Btot  and  the  typical  support 
of  each  plot,  given  by  e.g.  its  full  width  at  half  maximum. 

Now  substitute  (9.17)  into  (9.4)  to  give  the  range  resolution  of  the  chirped  pulse: 


A  r  ~ 
— res  — 


C 

2  (Hr  +  1/r)  ' 


(9.18) 


Given  a  pulse  width  r  and  a  desired  range  resolution  Arres,  we  can  find  the  necessary 
chirp  parameter  /x  from  (9.18): 


i)i 


(9.19) 


There  are  two  values  possible  here,  corresponding  to  an  “up  chirp”  ( \i  >  0)  and  a  “down 
chirp”  (fj.  <  0). 

The  chirped  pulse  has  an  effective  pulse  width  of 

reff  =  T> —  =  ~T~\ — T  i  /  =  |  |  2  i  i  '  (9.20) 

mr  +  l/T  \v\t2  +  1 

Alternatively,  its  compression  ratio  is  r/reff  =  |/x|  r2  +  1. 

Choosing  a  Sampling  Interval  for  a  Chirped  Pulse 

We  require  our  pulse  sampling  procedure  to  capture  the  emitted  wave’s  modulation  ade¬ 
quately,  but  there  is  no  point  in  making  the  sampling  frequency  excessively  high.  Given 
that  we  will  sample  at  a  constant  rate,  then  when  the  emitted  signal  is  chirped,  we  need  to 
sample  at  two  or  more  times  the  largest  modulation  frequency  present  to  satisfy  Nyquist’s 
criterion  for  detecting  this  modulation.  However,  if  we  do  this  and  then  reduce  the  chirp 
to  zero  (so  that  its  largest  modulation  frequency  goes  to  zero),  the  sampling  frequency 
will  then  go  to  zero:  we  will  stop  sampling  it!  To  prevent  this,  we  also  require  to  sample 
the  pulse  at  least,  say,  10  times.  So  we  must  set  a  lower  bound  on  the  sampling  frequency. 

9  Convolving  two  sequences  is  exactly  the  same  procedure  as  setting  the  reverse  of  one  sequence  to  be 
a  moving-mean  template  of  weights  that  “smoothens”  or  “smudges”  the  other  sequence.  Which  sequence 
plays  which  role  here  is  immaterial,  because  A  *  B  =  B  *  A. 
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To  quantify  this  further,  note  that  (9.5)  gives  the  chirp  modulation  as  el7r/ii2,  so 
the  instantaneous  non-angular  frequency  of  this  modulation  at  time  t  is  approximately, 
from  (7.6), 

L^t2  =  Mt.  (9.21) 

Irrespective  of  whether  /i  is  positive  or  negative,  the  absolute  value  of  the  modulation 
frequency  is  \n\t,  which  attains  a  maximum  value  of  \n\ r  at  the  end  of  the  chirp.  We 
sample  at  7q  times  this  rate  at  the  very  least  (where  rq  is  at  least  the  2  of  the  previous 
paragraph  to  satisfy  Nyquist);  i.e.  the  sampling  frequency  must  be  at  least  rq|/t|r.  But 
also,  we  require  at  least  n2  samples  in  the  whole  pulse  (where  n2  was  10  in  the  previous 
paragraph),  giving  a  sampling  frequency  of  at  least  n2/r.  So  choose  the  larger  of  7q|/t|r 
and  n2/r;  any  larger  number  would  only  lead  to  wasteful  oversampling.  Now,  sampling 
frequencies  are  reciprocals  of  sampling  times,  so  this  choice  of  larger  sampling  frequency 
corresponds  to  a  choice  of  smaller  sampling  time: 

ts  =  mini  — i— ,  —1 
t?q|/i|T  n2J 


(9.22) 


9.2  Range  Resolution  of  a  Barker-coded  Pulse 


The  bandwidth  of  a  binary  phase-coded  pulse  such  as  a  Barker  approximately  equals  the 
bandwidth  of  one  of  its  chips.10  Since  each  chip  is  just  the  carrier  modulated  by  a  top  hat, 
the  chip’s  bandwidth  equals  the  reciprocal  of  its  duration.  If  the  pulse  width  is  r,  each 
chip  of  an  n-chip  Barker  code  has  duration  r/n;  hence  the  bandwidth  of  the  pulse  must 
be  B  =  n/r,  giving  the  pulse  a  range  resolution  of 


A  r  ^ 
res 


C 

2B 


CT 

2  n 


(9.23) 


Compare  this  to  (9.3):  the  number  of  chips  n  plays  the  same  role  here  as  the  compression 
ratio  of  a  chirped  pulse.  Alternatively,  the  Barker’s  effective  pulse  width  is 


TeS  =  l/B  =  T/n, 

giving  the  Barker- n  pulse  a  “compression  ratio”  of  t/tsq  =  n. 


(9.24) 


10  Representative  Range— Doppler  Plots 

In  this  section  we  generate  range-Doppler  plots  for  various  coherent  processing  intervals 
of  non-chirped  and  chirped  pulses.  The  first  plots  use  32  pulses  per  interval,  first  without 
and  then  with  Doppler  windowing.  The  number  of  pulses  per  CPI  in  the  windowed  set  is 
then  increased  to  256  as  a  simple  example  of  showing  how  this  increase  leads  to  a  better 
Doppler  resolution — but  not  a  better  range  resolution.  All  plots  use  a  10  GHz  carrier. 

10This  statement  can  be  verified  numerically  when  the  chip  polarities  are  chosen  randomly,  but  I  doubt 
that  any  solid  proof  of  it  is  possible,  because  bandwidth  itself  is  not  precisely  defined.  Clearly  it  will  fail 
if  the  chips  all  have  the  same  polarity.  Its  truth  is  probably  due  to  the  idea  that  generally  bandwidth 
is  inversely  proportional  to  signal  length,  so  that  a  short-duration  feature  of  the  signal  (such  as  a  chip) 
requires  more  bandwidth  than  a  long-duration  feature,  and  so  the  bandwidth  of  that  short-duration  feature 
will  account  for  most  of  the  total  bandwidth. 
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10.1  Series  of  Non-Chirped  Rectangular  Pulses,  Before  and 
After  Windowing 

For  a  target  at  500  m  approaching  at  100  m/s,  the  information  on  selecting  valid  pulse 
parameters  in  Section  8  allows  us  to  set  parameters  for  a  series  of  rectangular  non- 
chirped  pulses.  Equation  (8.8)  specifies  that  the  PRI  should  lie  between  3.3  and  75  mi¬ 
croseconds,  so  we  have  used  the  average  value  of  39  microseconds,  corresponding  to  a 
PRF  of  about  26  kHz.  A  range-Doppler  plot  with  no  windowing  is  shown  at  the  left  of 
Figure  14.  The  peak  correlation  occurs  at  range  510  metres  and  velocity  —105  m/s. 


Target  distance: 

500  m 

Target  recession  velocity: 

—  100  m/s 

Range-bin  width: 

Arres  =  150  m 

Pulse-train  length: 

N  =  32  pulses  per  CPI 

Pulse  width: 

T  =  1  flS 

PRI  (8.8): 

T  =  39  gs  (PRF  =  26  kHz) 

Blind  range  and  max.  unambiguous  range  (8.1): 

150  m  and  5846  m  respectively 

Max.-detectable  speed  and  resolution  (8.10),  (8.11): 

192  m/s,  12  m/s 

Chirp  factor  (9.19): 

g  =  0  MHz2  (reff  =  1  gs) 

Sampling  interval  (9.22): 

ts  =  0.067  gs  (n-x  =  2,  n2  =  15) 

Emitted  signal  (modulation) 

Received  signal  (modulation) 

0. 

0. 

0. 

tT  3  n 

os  “■ 

^  o. 

o  ° 

if  o- 

o. 

o. 
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Figure  13:  A  CPI  of  32  pulses  with  range  resolution  of  10  metres,  no  chirp, 

PRI  =  39  gs  (PRF  =  26  kHz),  pulse  width  1  gs,  sampling  interval  0.067  /is, 
and  no  Doppler  windowing.  Left:  modulation  of  emitted  signal.  Right:  re¬ 
ceived  signal. 

At  the  bottom  of  Figure  14  is  a  range-Doppler  plot  made  of  the  same  target,  but 
now  incorporating  Doppler  windowing.  The  Doppler  tail  is  now  completely  absent  at  the 
expense  of  a  slight  broadening  of  the  Doppler  peak.  The  peak  correlation  now  occurs  at 
range  500  metres.  The  velocity  still  peaks  at  —105  m/s. 

This  “tidying  up”  of  the  range-Doppler  plot  is  the  reason  why  windowing  functions  are 
useful,  especially  when  windowing  is  used  in  the  presence  of  background  noise  and  other 
targets. 
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|  Frequency  spectrum  |,  32  pulses,  p  =  0  MHz2 
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Figure  1 f:  Top:  The  CPI  of  32  pulses  in  Figure  13  gives  this  range- 
Doppler  plot.  No  Doppler  windowing  has  been  used.  Bottom:  The  same 
scenario  with  Doppler  windowing  included. 
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|  Frequency  spectrum  |,  256  pulses,  p  =  0  MHz2 


Figure  15:  A  CPI  with  the  same  parameters  used  for  Figure  If,  with 
Doppler  windowing,  but  now  using  256  pulses.  The  higher  number  of  pulses 
increases  the  Doppler  resolution. 


Figure  15  shows  range-Doppler  for  the  same  parameters  used  for  Figure  14,  including 
windowing,  but  now  with  256  pulses  per  CPI.  The  range  result  is  again  500  metres,  but 
the  higher  number  of  pulses  has  increased  the  Doppler  resolution,  so  that  the  velocity  is 
now  much  more  localised  near  its  peak  at  —102  m/s. 
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10.2  Series  of  Chirped  Rectangular  Pulses  with  Doppler 
Windowing 

The  parameters  in  the  table  below  have  been  used  to  create  a  series  of  chirped  pulses, 
again  for  a  target  at  500  m  approaching  at  100  m/s.  Equation  (8.8)  specifies  that  the 
PRI  should  lie  between  3.3  and  75  microseconds,  so  we  have  used  the  average  value  of 
39  microseconds. 


Target  distance: 

Target  recession  velocity: 

Range-bin  width: 

Pulse-train  length: 

Pulse  width: 

PRI  (8.8): 

Blind  range  and  max.  unarnb.  range: 
Max.-detectable  speed  and  resolution  (8.10),  (8.11): 
Chirp  factor  (9.19): 

Sampling  interval  (9.22): 


500  m 
—  100  m/s 
Arres  =  10  m 
N  =  32  pulses  per  CPI 

T  =  1  fiS 

T  =  39  [is  (PRF  =  26  kHz) 

150  m  and  5846  m  respectively 
192  m/s,  12  m/s 
/r  =  14.0  MHz2  (reff  =  1/15  //s) 
ts  =  0.036  /rs  (rii  =  2,  n2  =  15) 


The  effective  pulse  width  is  reg  =  1/15 /us,  giving  a  chirp  factor  of  /i  =  14  MHz2.  Some 
intuitive  feel  for  the  value  of  /r  results  from  writing  it  as  14  MHz/gs.  Equation  (9.6)  gives 
the  carrier  frequency  of  the  emitted  pulse  as  /  =  /o  +  /A,  so  for  the  1  /j, s  duration  of  the 
chirp  the  carrier  frequency  increases  at  a  rate  of  [i  =  14  MHz/^rs,  giving  a  total  frequency 
increase  of  14  MHz.  This  is  the  swept  bandwidth  of  one  pulse.  (This,  of  course,  is  just  a  re¬ 
statement  of  (9.16),  which  gives  the  frequency  change  over  a  pulse  width  as  Bswept  =  |//| r.) 
Compare  this  14  MHz  frequency  increase  with  the  carrier  frequency  of  10  GHz:  the  signal 
is  indeed  narrow  band. 

Figure  16  (top)  shows  the  resulting  range-Doppler  plot.  The  peak  correlation  occurs  at 
range  of  498  metres  and  a  velocity  of  —105  m/s.  The  chirp  gives  a  better  range  resolution 
than  the  non-chirped  pulse  trains  in  Figures  14  and  15,  although  the  fact  that  we  don’t 
quite  get  a  peak  at  500  metres  is  probably  just  due  to  a  chance  placement  of  the  sampling. 
The  bottom  plot  of  Figure  16  shows  the  same  scenario  but  now  using  256  pulses,  which 
gives  better  Doppler  resolution. 
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|  Frequency  spectrum  |,  32  pulses,  p  =  14  MHz2 


|  Frequency  spectrum  [,  256  pulses,  p  =  14  MHz2 


Figure  16:  Top:  A  CPI  of  32  pulses  with  Doppler  windowing,  range  res¬ 
olution  of  10  metres,  chirp  p  =14-0  MHz2,  PRI  =  39  ps,  pulse  width  1  ps, 
sampling  interval  0.036  ps.  The  chirp  increases  the  range  resolution.  Bot¬ 
tom:  The  same  for  a  CPI  of  256  pulses,  showing  a  clear  increase  in  Doppler 
resolution. 
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11  Creating  Range  and  Velocity  Offsets 
for  Use  in  Jamming 

A  radar  target  might  actively  record  the  signal  and  re-emit  a  copy,  perhaps  altered,  after 
some  time  interval.  The  alteration  is  accomplished  by  a  set  of  electronics  called  a  tap, 
which  increases  the  signal  strength  (wave  amplitude),  and  possibly  introduces  a  time- 
dependent  phase  shift.  We  describe  these  alterations  in  this  section.  As  in  Section  2.3,  the 
analytic  signal  emitted  by  the  radar  at  time  t  is  u(t)  eluJ()t,  and  the  target  motion  is  also 
modelled  as  in  that  section.  The  question:  how  should  the  tap  alter  the  signal  to  offset 
its  range  and  velocity  on  the  receiver’s  range-Doppler  map? 


11.1  Creating  a  Range  Offset 


The  tap  delay  Ttap  that  gives  a  desired  range  offset  is  straightforward  to  compute.  Delaying 
a  tap’s  return  by  Ttap  leads  the  radar  to  measure  the  tap  to  lie  a  distance  of  A r  =  cTtap/2 
further  away  than  its  true  value.  So  for  any  given  A r,  the  required  delay  is 


^taP  -  2  A r/c . 


(11.1) 


For  example,  if  we  wish  the  target  to  appear  to  be  150  metres  further  away  than  what  it 
really  is,  set  A r  =  150  m  in  (11.1)  and  infer  that  a  tap  must  delay  the  signal  for  a  time  of 


2  x  150  m 
tap  300  m//xs 


(11.2) 


11.2  Creating  a  Velocity  Offset 

A  tap  can  introduce  a  velocity-offset  structure  to  the  range-Doppler  plot  by  injecting  a 
time-dependent  phase  change  to  each  pulse  it  returns.  We  can  see  this  in  the  following 
way.  Refer  to  the  third  option  of  (2.22): 

gc(t)  ~  A{rG)u{t-td)e^e-^-c^+^).  (11.3) 

Suppose  that  a  tap  multiplies  each  I/Q  number  it  generates  by  e~^o2/cAvt ^  which  has  the 
effect  of  replacing  v  in  (11.3)  with  v  +  Av.  We  then  expect  a  new  velocity  of  v  +  Av  to 
be  measured.  In  practice,  of  course,  the  tap  has  a  clock  of  its  own  which  is  offset  from  the 
receiver’s  clock  by  some  toffset.  So  replace  the  factor  of  e~uoo2!c Avt  by 

e~luJ°c  AU*  +  £offset)  —  e-M)c  Avt  x  a  constant  phase  factor.  (11.4) 


So  this  tap’s  clock  offset  produces  a  constant  phase  shift,  which  is  of  no  real  consequence 
in  the  radar  signal  processing.  The  additional  Doppler  shift  detected  by  the  radar  that 
arises  from  Av  follows  from  (2.26),  whose  linearity  implies  that 


A  ujd  =  —u02A  v/c. 


(11.5) 


Here  we  have  a  recipe  for  creating  a  single  Doppler  offset  from  a  given  tap.  But  in  addition, 
although  each  tap  is  assumed  capable  of  returning  only  one  time-delayed  signal  to  the  radar 
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for  each  pulse  incident  on  it—  producing  a  single  major  range  offset — the  tap  is  assumed 
to  be  capable  of  building  this  return  by  summing  multiple  altered  copies  of  the  incident 
pulse,  where  each  alteration  consists  of  an  amplitude/phase  variation  made  to  the  pulse. 
Each  of  these  copies  that  are  returned  simultaneously  will  give  rise  to  its  own  peak  in  the 
Doppler  spectrum  produced  by  the  radar,  with  all  of  these  peaks  at  the  same  range  for  a 
given  tap. 

Returning,  say,  two  altered  copies  of  the  incident  signal  back  to  the  radar  from  tap  n, 
corresponding  to  velocity  offsets  Aiq  and  Av2,  effectively  requires  multiplying  one  copy 

by  Wn^  e~luJ°2/cAvit  [  =  Wn^  elAuD’nt  ]  and  another  copy  by  Wn"1  e~luJ°2^c/^V2t  [  =  Wn^  exAu°<nt 

for  weights  Wn  \  Wn  ^ ,  and  summing  these  two  copies.  This  corresponds  to  multiplying  the 
incident  signal  by  the  single  number 

e~iuJ°2/cAvit  +  e-^oVcAv^t  j-  =  wWeiAJ^nt  +  w(2) eiAu^„t j  _  (n.6) 

This  new  factor  becomes  a  generalised  “Doppler  offset”  factor  in  the  I/Q  returned  by  the 
tap.  It  multiplies  the  right-hand  sides  of  the  expressions  in  (2.34)  and  (2.35).  It  needn’t 
have  unit  magnitude.  It  will  produce  two  peaks  in  the  receiver’s  range-Doppler  map  at 
different  Doppler  values  but  at  a  single  range  value.  Of  course,  in  general  there  can  be  any 
number  of  terms  in  the  phase  sum.  In  fact  we  can  synthesise  an  arbitrary  Doppler  profile 

(?)  iAo/-^  t 

for  tap  n  by  summing  some  possibly  large  number  of  Wn  'e  D>n  terms,  although  whether 
the  tap  is  capable  of  generating  this  large  number  of  terms  depends  on  its  hardware. 

In  practice,  only  absolute  values  of  the  Doppler  spectral  components  are  conventionally 

plotted  in  each  range  bin  to  build  the  range-Doppler  map;  phase  information  is  discarded. 

(i) 

But  inserting  arbitrary  time-independent  phases  into  the  weights  w)i  should  make  no 
difference  to  the  synthesised  range-Doppler  plot,  and  indeed  they  do  not,  as  noted  af¬ 
ter  (11.10)  ahead. 

So  if  we  require  to  synthesise  the  Doppler  profile  in  range  bin  “n”  of  an  already  existing 
range-Doppler  plot  that  was  formed  using,  say,  32  pulses/CPI  (and  so  has  32  velocity 
values  plotted  at  each  range),  then  we  might  be  able  to  use  32  velocity  offsets:  the  values 

of  A ojd  n  to  A ujyD  ^  will  be  those  values  corresponding  to  the  velocities  on  the  range- 

(1)  (32) 

Doppler  plot  by  way  of  (11.5),  and  their  weights  Wn  to  Wn  ;  will  be  proportional  to  the 
values  being  plotted  on  the  “z  axis”  of  the  range-Doppler  plot.  But  given  that  there  is 
always  some  velocity  broadening  along  the  Doppler  axis,  we  might  well  need  to  specify 
more  than  32  velocity  offsets  to  create  the  desired  Doppler  profile. 


The  Taps’  Impulse  Response 


The  ship’s  impulse  response  appeared  in  the  calculation  of  Section  3.1.  Here  we  show  how 
to  modify  that  calculation  to  model  the  taps’  impulse  response  for  including  range  and 
velocity  offsets.  In  particular,  the  return  is  delayed  as  per  (11.1)  (now  using  a  subscript  n 
for  tap  n),  so  that  a  delayed  return  from  tap  n  increments  element  number  kn  of  the 
impulse-response  array  h,  where  kn  was  calculated  using  (3.7)  for  no  delay,  but  is  now 
calculated  from 


kn  —  1  =  round 


n  T  A  Tn  -Rnear) 

Cts 


(11.7) 
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|  Frequency  spectrum  |,  32  pulses,  p  =  9  MHz2 
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Figure  17:  A  two-tap  synthesis  of  two  range  bins,  using  the  parameters 
in  (11.9)  and  (11.10).  The  first  range  bin  has  four  velocity  offsets  and  the 
second  has  two  velocity  offsets. 


The  return  itself  now  incorporates  the  weighted  sum  of  Doppler  offsets,  such  as  in  (11.6). 
For  tap  n,  (3.8)  becomes 


modulation  of  return 


^  ,  jt(^t(Tn  e^O.n  e~iaJ0§  x  ran§e  at  emission  ^  W<J)  e~iuJ0  §  Avjt 

4?r r2  y  £0c  . 

^  j^tGtO'n'  gi tj>0n  e-iw0  §  X  range  at  emission  wfj) 

4-71-7-2  y  £qC 

(11.8) 


where  t  is  the  time  at  the  start  of  the  current  pulse  being  emitted. 


11.3  Example  of  Synthesising  Range  and  Doppler  Offsets 

Figure  17  shows  an  example  of  a  range-Doppler  profile  built  from  a  set  of  velocity  off¬ 
sets  A  that  give  rise  to  corresponding  Aw^n,  using  two  taps,  with  Blackman  windowing 
in  Doppler,  and  incorporating  an  arbitrarily  chosen  signal-to-noise  ratio  of  20  dB,  as  shown 
in  the  following  table.  A  10  GHz  carrier  is  used. 


Target  distance: 

Target  recession  velocity: 

Range-bin  width: 

Pulse-train  length: 

Pulse  width: 

PRI  (8.8): 

Blind  range  and  max.  unarnb.  range: 
Max.-detectable  speed  and  resolution  (8. 
Chirp  factor  (9.19): 

Sampling  interval  (9.22): 


500  m 
0  m/s 

Arres  =  15  m 
N  =  32  pulses  per  CPI 

T  =  1  pS 

T  =  40  ps  (PRF  =  25  kHz) 
150  m  and  5996  m  respectively 
,(8.11):  187  m/s,  12  m/s 

p  =  9  MHz2  (reff  =  1/10  ps) 
ts  =  0.1  ps  (ni  =  2 ,  n2  =  10) 
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Tap  1  offsets  range  by  Ar1  =  —340  m  from  the  nominal  target  range  of  500  m,  and 
has  four  Doppler  offsets,  at  velocities  of  135,  88,  50,  and  —24  m/s  from  the  nominal  target 
recession  velocity  of  0  m/s,  using  tap  weights  of  1,  0.5,  0.2,  0.8  respectively: 


The  parameters  for  tap  2  are 


Ar 

L  =  —340  m, 

Av 

f ,  Auf>] 

>  =  {135,  88 

,  50,  -24} 

'  (1)  (2) 
w{  w\ 

(3)  (4) 1 

>  =  {L  0.5, 

0.2,  0.8}. 

2  are 

II 

<1 

—  100  m, 

{a41}, 

> 

to 

II 

{-50,  100} 

m/s, 

42)}  = 

{1,  2}- 

(11.9) 


(11.10) 


As  expected,  if  we  attach  a  random  time-independent  phase  to  each  of  the  tap  weights, 
these  plots  are  visually  unaffected. 


12  Final  Comments  and  Acknowledgements 

This  report  has  addressed  two  aims.  The  first  has  been  to  act  as  a  tutorial  in  the  funda¬ 
mentals  of  radar  signal  processing.  My  approach  has  been  to  start  from  the  first  principles 
of  electromagnetic  theory,  build  an  expression  for  the  I/Q  signal  that  is  returned  by  an 
arbitrary  number  of  scatterers,  then  show  how  this  yields  range  and  Doppler  information 
when  correlated  with  the  emitted  signal.  I  have  also  laid  out  some  of  the  theoretical  rea¬ 
sons  behind  the  various  bounds  placed  on  standard  waveform  parameters  such  as  pulse 
width  and  chirp  bandwidth. 

In  the  appendices  that  follow,  I  have  analysed  the  Hilbert  transform  to  show  where 
it  comes  from  and  why  it’s  useful.  My  approach  stresses  how  it  can  be  constructed  in  a 
proper  mathematical  way  without  the  arbitrariness  found  in  the  standard  literature,  where 
the  transform  is  typically  defined  as  a  principal  value  for  no  apparent  reason  other  than 
that  this  avoids  a  divergent  integral.  My  approach  does  produce  this  same  definition,  but 
now  the  principal  value  arises  quite  naturally. 

I  have  also  explored  a  pen-and-paper  approach  to  convolution  and  correlation  in  an 
appendix.  The  fact  that  “long  multiplication”  is  actually  a  convolution  is  well  known,  but 
I  hope  the  appendix  gives  additional  insight,  in  that  it  employs  a  little-known  form  of 
long  multiplication  that  matches  the  convolution  routine  far  more  closely  than  the  usual 
“school”  method  of  long  multiplication  does. 

The  second  aim  of  this  report  has  been  to  describe  the  mathematics  of  how  a  signal 
can  be  modified  by  a  jammer  so  as  to  produce  false  targets  in  the  radar  receiver. 

I  wish  to  thank  Roland  Keir  for  his  close  involvement  in  the  writing  of  this  report. 
Its  contents  also  benefitted  from  feedback  and  discussions  with  Gavin  Dickeson,  Len  Hall, 
Stephen  Howard,  and  Warren  Marwood. 
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Appendix  A  Some  Discussion  of  the  Hilbert 

Transform 


This  appendix  gives  a  derivation  of  the  Hilbert  transform  with  substantially  more  mathe¬ 
matical  detail  than  typically  found  in  radar  texts.  It  is  mostly  my  own  analysis  of  why  the 
Hilbert  transform  takes  the  form  that  it  does.  I  have  yet  to  find  a  radar  textbook  that  dis¬ 
cusses  the  Hilbert  transform  in  a  way  that  as  a  mathematical  physicist  I  find  illuminating; 
hence  I  have  taken  an  approach  here  that  does  speak  to  the  mathematical  physicist.  But 
additionally,  the  analysis  here  explains  why  the  standard  approach  of  radar  texts  works, 
when  that  approach  “fixes”  a  divergent  integral  in  an  apparently  ad  hoc  manner. 

Given  a  real  signal  s(t),  we  are  generally  not  interested  in  its  carrier;  the  information 
in  the  signal  is  contained  in  its  running  amplitude  and  phase.  The  amplitude  and  phase 
conventionally  define  the  phasor  representation  of  the  signal.  We  will  set  s(t)  to  be  the 
real  part  of  a  complex  signal  sc(t),  and  then  determine  the  choice  of  imaginary  part  that 
gives  a  phasor  clearly  embodying  the  signal’s  amplitude  and  phase.  The  resulting  special 
choice  of  complex  signal  will  be  called  the  analytic  signal. 

Suppose,  first,  that  our  real  signal  is  s(t)  =  cos  5 1.  Set  this  to  be  the  real  part  of  sc(t); 
we  desire  to  have  sc(t)  =  et54,  since  this  represents  a  phasor  of  unit  length  (corresponding 
to  the  unit  amplitude  of  cos  5 1),  and  constant  spin  rate  u  =  5,  corresponding  to  the  an¬ 
gular  frequency  co  =  5  of  the  real  signal.  An  alternative  choice  is  sc(t)  =  e”*54,  which  also 
has  real  part  cos  5 1.  By  convention,  we’ll  set  sc(t)  to  have  positive  frequency:  sc(t)  =  el54. 
Notice  that  this  contains  one  “phasor”  frequency  co  =  5  with  unit  weight.  Compare  this  to 
the  real  signal  s(t)  =  cos  5 1  =  !/2ei54  +  i^e-154,  which  is  built  from  two  Fourier  frequen¬ 
cies  co  =  ±5,  each  with  weight  1/2. 

This  construction  of  a  useful  complex  signal  from  the  real  one  is  trivial  in  this  case, 
but  how  might  we  go  about  it  given  a  more  complicated  s(t)?  After  all,  we  can  always 
write  s(f)  =  cos  5 1  in  arbitrary  ways,  such  as 

s(t)  =  A(t )  cos(3 1  +  sin  t)  =  B(t)  cos(12f  +  e4)  (Al) 


for  some  real  A(t)  and  B(t),  but  that  does  not  imply  that  it  will  be  useful  to  define 

sc(t)  =  A(t)  exp  i(3f  +  sin  t) ,  or  sc(t)  =  B(t )  exp  i(12 1  +  e4) .  (A2) 

That  these  two  choices  of  sc(t)  are  generally  not  equal  can  be  shown  as  follows.  The  real 
parts  of  the  two  choices  are  equal  by  construction,  so  we  must  show  that  the  choices’ 
imaginary  parts  are  generally  not  equal.  Do  this  by  using  (Al)  to  solve  for  A(t)  and  B(t): 


A(t)  = 


s{t) 


,  B(t)  = 


s(t) 


(A3) 


cos(3t  +  sinf)  ’  w  cos(12t  +  e4) 

The  imaginary  parts  of  the  two  choices  of  sc(t)  in  (A2)  must  then  be,  from  (A3), 

A(t)  sin(3f  +  sin  t)  =  s(t)  tan(3 1  +  sin  t) ,  (A4) 

and 


B(t)  sin(12t  +  e4)  =  s(t)  tan(12 1  +  e4) .  (A5) 

These  two  imaginary  parts  are  generally  not  equal,  so  the  two  choices  of  sc(t)  in  (A2)  are 
also  generally  not  equal.  And  because  the  functions  3t  +  sinf  and  12t  +  e4  were  arbitrarily 
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chosen,  we  see  that  there  are  an  infinite  number  of  choices  of  sc(t)  that  all  have  the  same 
real  part  s(t),  but  which  have  different  imaginary  parts.  The  question  is,  which  choice 
of  sc(t)  might  be  “best”? 

Different  choices  of  sc(t)  will  generally  have  frequency  spectra  that  don’t  resemble 
the  spectrum  of  s(t)  =  cos  5 t.  We  define  the  best  choice  of  sc(t)  for  s(t)  =  cos  5 1  to  be 
the  standard  phasor:  sc(t)  =  e1"5*.  The  frequency  spectrum  of  this  sc(t)  has  no  nega¬ 
tive  frequencies,  and  has  twice  the  amount  of  the  positive  frequency  that  appears  in 
s(t)  =  cos  5 t  =  1/2  el5t  +  This  idea  can  now  be  applied  to  any  real  signal  s(t): 

choose  sc(t)  to  be  the  function  whose  frequency  spectrum  has  no  negative  frequencies, 
and  whose  positive  frequencies  are  the  same  as  those  of  the  real  signal  s(t),  but  with 
twice  the  weighting.  Such  a  function  is  easily  constructed:  Fourier-transform  s(t)  to  give 
its  spectrum,  remove  the  negative  part,  double  the  positive  part,  and  take  the  inverse 
Fourier  transform,  resulting  in  sc(t).  This  special  choice  of  sc(t)  is  called  the  analytic 
signal  corresponding  to  the  real  signal  s(t).  In  Fourier  language,  if  the  real  signal  s(t)  has 
spectrum  S(uj): 

/OO 

S(ui)  elujt  dcu  ,  (A6) 

-OO 

then  the  analytic  signal  is 

r  00 

8c(t)  =  /  2 S(u)  dw .  (A7) 

Jo 

Before  evaluating  this  integral,  we’ll  give  an  alternative  discussion  of  the  analytic  signal 
that  also  leads  to  (A7). 

A.l  Alternative  Derivation  of  the  Analytic  Signal 

Here  is  another  approach  to  constructing  a  “best”  choice  of  complex  signal  sc(t)  that  is 
not  related  to  the  discussion  above.  Given  a  real  signal  s(t),  suppose  we  construct  a  linear 
operator  C  that  converts  s(t)  to  a  complex  signal  by  adding  an  imaginary  part  i s(t),  where 
s(t)  is  some  special  real  function  to  be  determined: 

C{s(t)}  =  s(t)  +  is(t) .  (A8) 

We  require  only  that  for  positive  frequencies  H  (we  are  taking  a  minimalist  approach  by 
dealing  with  positive  frequencies  only),  C  converts  cos  Qt  and  sin  ilt  into  phasors,  viz.,  the 
complex  exponential  e1  ,  possibly  with  some  scale  factor.  Specifically, 

CjcosUt)  =  elilt ,  (A9) 

C{sin  fit}  =  aelQt  for  some  a  .  (A10) 

These  three  requirements  (A8)-(A10)  will  suffice  to  determine  C  fully.  Write  a  =  a  +  ib 
for  a,  b  real,  and  apply  (A8)  to  the  left-hand  side  of  (A10)  to  give 

sin  f it  +  isin  Qt  =  (a  +  ib)  (cos  Df  +  isin  f It) .  (AH) 

Equating  the  real  and  imaginary  parts  of  (All)  and  making  use  of  the  linear  independence 
of  sine  and  cosine  leads  to  a  =  — i,  in  which  case  sin  fit  =  —  cos  Ilf.  Thus 

Clsinflf}  =  —  ieint .  (A12) 
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We  now  have  enough  information  to  determine  how  C  acts  on  sinusoids  of  both  positive 
and  negative  frequencies  uj.  Use  the  fact  that  cosine  is  an  even  function  (cos#  =  cos  |#|) 
and  sine  is  odd  (sin#  =  sgn#  sin  |#|)  to  write 


C{cosu;t}  =  C{cosMt}  =  elM\ 

C{sinu;t}  =  Cjsgna;  sin  |w|t}  —  isgncu  , 

and  therefore,  on  adding  the  first  equation  above  to  i  x  the  second  one, 

C{eiujt}  =  (1  +  sgncu)  e1Mt  =  2#(w)  eiujt , 
where  0(uj)  is  the  unit  step  function:11 

!1 ,  uj  >  0 
1/2  ,  uj  =  0  . 

0  ,  uj  <  0 


(A13) 


(A14) 


(A15) 


We  can  now  use  the  linearity  of  C  to  build  C'js(t)}  for  a  general  real  signal  s(t).  Begin 
with  the  Fourier  decomposition  (A6),  writing 

/oo  roo  r  oo 

S(uj)C{eiujt}duj  =  /  S(cj)20(cj)eiujtdu=  /  2S(uj)  eiut  du .  (A16) 

-oo  J— oo  J  0 

But  the  last  integral  above  is  just  the  right-hand  side  of  (A7),  which  shows  that  C'{s(t)}  is 
identical  to  the  analytic  signal  sc(t)  of  (A7).  So  this  alternative  approach  to  complexi¬ 
fying  s(t)  has  reproduced  the  analytic  signal  described  in  (A1)-(A7).  We  need  now  only 
evaluate  (A7). 

Notice  that 

__ — - — _  (A9) 

cos  f it  =  sin  f It  =  cosfflt  —  7r/2) ,  and 
_ _ _  (A.12) 

sin  fit  =  — cos  fit  =  sinfflt  —  7r/2) .  (AIT) 

This  means  that  the  imaginary  part  of  the  complex  signal  built  from  cos  fit  is  found  by  simply 
delaying  cos  fit  by  90°.  The  same  is  true  for  the  imaginary  part  of  the  complex  signal  built 
from  sin  fit.  Because  Fourier  analysis  states  that  a  signal  can  be  considered  as  a  sum  of  sines  and 
cosines,  we  conclude  that  the  complex  part  of  the  analytic  signal  is  built  by  taking  each  of  the 
sinusoidal  components  of  the  real  signal  and  delaying  them  by  90°,  or  a  quarter  period. 


A. 2  Evaluating  the  One-Sided  Fourier  Integral  (A7) 

To  evaluate  (A7),  begin  by  incorporating  Fourier  theory  into  the  current  analysis  through 
writing  (A7)  as  an  integral  from  —  oo  to  oo: 

sc(t)=  f°°2S(uj)  eiujt  dw  =  f°°  9{uj)  2S{uj)  elwt  du;  =  H  S(uj)  eiwt(l  +  sgncu)  do; 

JO  J— oo  J —oo 

=^=^=  s(t)  +  f  S (uj)  elujt  sgn  uj  duj .  (A18) 

J —oo 

nI  use  0(0)  =  1/2,  corresponding  to  the  choice  sgnO  =  0.  This  aligns  with  the  sense  of  continuity  of 
Fourier  sums  within  the  Fourier  transform  theory  of  discontinuous  functions. 
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The  sgncu  in  the  last  integral  of  (A18)  can  be  converted  to  another  Fourier  integral  through 
the  following  piece  of  gymnastics.  First,  use  the  fact  that 


sm  ojt 

- dr 

r 


7 r  sgn  uj  . 


(A19) 


Although  sm^T  ig  not  defined  for  r  =  0,  it  has  a  “removable  discontinuity”  there:  a  single 
“hole”  in  the  curve  depicting  the  function,  as  it  were.  You  will  sometimes  see  it  written 
that  =  1  when  x  =  0.  But  this  is  not  so;  rather,  only  the  limit  exists:  limx^o  =  1, 
and  distinguishing  between  these  two  statements  is  precisely  why  the  concept  of  a  limit  was 
invented  three  centuries  ago.  The  term  “removable  discontinuity”  describes  the  fact  that 
functions  with  isolated  holes  in  their  otherwise  continuous  curves  still  have  well-defined 
integrals.  So  the  presence  of  a  single  “hole”  in  the  curve  of  the  sine  function  doesn’t  stop 
the  integral  (A19)  from  being  well  defined.  But  that  means  its  principal  value  is  also  well 
defined: 


f°°  sin  cur 

I-  oo  t 


dr  =  7r  sgn  uj  , 


(A20) 


where  the  principal  value  is  defined  as 


lim 

£->0+ 


/(x) dx  + 


f(x) dx 


(A21) 


(For  example,  whereas  f^dx/x  is  undefined,  the  principal  value  -j-^^dx/x  equals  zero 
because  1/x  is  an  odd  function.) 

Writing  (A19)  in  principal-value  form  (A20)  allows  us  to  introduce  the  following 
principal-value  integral: 


f°°  cos  cot  . 

+  - dr  =  0 

J — oo  T 

(for  which  the  principal  value  is  required )  to  conclude  that 


(A22) 


f°°  glWT  /■  o°  COS  LOT  +  i  SU1  COT 

+  - dr  =  +  - dr  =  itt  sgn  uj  . 

J— oo  "7”  ./—oo  T 


(A23) 


Equation  (A23)  is  almost  a  Fourier  integral — but  only  in  the  principal-value  sense.  We 
can  eliminate  the  need  for  a  principal  value  and  convert  (A23)  fully  to  a  Fourier  integral 
by  defining  the  following  continuous  generalised  function: 


cr(f) 


1/t  t  7^  0 
0  t  =  0. 


(A24) 


Defining  a(t)  as  continuous  makes  it  a  generalised  function  in  the  same  class  as  the  Dirac  delta 
function,  in  that  both  of  these  functions  can  only  be  realised  in  practice  as  limiting  cases  of 
sequences  of  functions.  Two  such  sequences  are  [4] 


n5(t) 


lim 

e— >0+ 


t2  +  £2 


r{t) 


lim 

e-M3+ 


t 2  +  £2 


(A25) 


whose  similarity  reveals  a  close  relationship  between  iv8(t)  and  a(t),  which  we’ll  see  more  of  shortly. 


Introducing  a(t)  enables  us  now  to  write  the  principal-value  integral  (A23)  as  a  Fourier 
integral  which  is  not  a  principal  value,  and  which  therefore  can  be  brought  into  the  fold 
of  Fourier  theory: 


elUTa(r)  dr 


( =  i.7rsgnu;) . 


(A26) 
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(This  is  an  important  point;  the  fact  that  aft)  is  continuous  and  yet  tr(0)  =  0  is  what 
converts  the  principal- value  integral  to  a  normal  integral.)  Now  the  last  integral  in  (A18) 
can  be  written  as  (with  do;  and  dr  repositioned  for  readability) 


/oo  _  1  r  oo 

S(ui)  elut  sgno;  du  =  —  /  dw5(w)  elaJH7rsgnw 

-oo  J — oo 


1  P  OO  P  oo 

=  —  /  d  w5(u)e“f  dreia;V(r) 

i-'TT  J — oo  —  oo 


(A27) 


Swapping  the  order  of  integration  yields 


/OO  1  P  oo  /*oo 

5(w)  eiwt  sgn  u  do;  =  —  /  dr  <t(t)  /  da;  5(w)  e^(t+r) 

-oo  "L’TT  J— oo  J— oo 

(A6)  1  A00  /  \  / .  |  \ 

—  dr  <t(t)  s(i  +  r) 

ITT  J  —oo 


1  /°° 

=  —  /  dr  cr(— r)  s(t  —  t) 

i-TT  J  —oo 

i  r°° 

=  —  dr  <r(r)  s(t  —  r) 

7T  J —oo 


=  -a(t)  *  s(t) ; 

7 r 


(A28) 

(A29) 

(A30) 

(A31) 


here  (A29)  follows  from  (A28)  by  essentially  just  a  change  of  variables  r  — >•  — r  (alter¬ 
natively,  this  expresses  the  fact  that  the  signed  area  under  any  function  is  unchanged  if 
the  function  is  reversed  left-right);  also,  (A30)  follows  from  (A29)  because  a  is  an  odd 
function;  and  (A31)  is  a  convolution.  Finally,  equation  (A18)  becomes 


sc(t)  =  s(t)  +  icr(t)/Tr  *  s(t) , 


(A32) 


which  we  compare  with  (A8)  to  write  the  Hilbert  transform  of  the  real  signal  as 


9-({s{f)}  =  s(t )  =  cr(t)/TT  *  s(t) . 


(A33) 


On  a  side  note,  it’s  apparent  that  any  one-sided  Fourier  transform  can  be  written  as  a  sum  of  a 
two-sided  Fourier  transform  and  a  Hilbert  transform.  We  can  see  this  by  writing  (A8)  as 

sc(t)  =  [l  + 17/"]  s(t) ,  or  !/2  sc(t)  =  1/i  [l  + 17/]  s(t) ;  (A34) 

now  call  on  (A6)  and  (AT)  to  write  the  last  expression  above  as 


(A35) 


Radar  texts  probably  universally  write  the  generalised  function  cr(t)  defined  in  (A24)  as 
simply  1/t.  Although  1/t  differs  from  a(t)  only  at  a  single  point,  there  is  a  marked  differ¬ 
ence  in  using  1/t  versus  using  aft)  in  Fourier  analysis:  for  example,  whereas  ffff  aft)  cos  t  d t 
exists  (and  equals  zero),  dt  does  not  exist.  Radar  texts  typically  invoke  the  prin¬ 
cipal  value  to  “fix”  the  divergence  of  an  integral  such  as  d t,  but  I  don’t  believe 

that  mathematics  should  be  done  by  ad  hoc  fixes  when  something  breaks;  there  is  no 
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a  priori  reason  for  why  enforcing  a  symmetrical  approach  to  a  singularity  in  an  integration 
should  be  a  reasonable  thing  to  do.  That  is,  one  should  not  invoke  a  principal  value  to 
“fix”  a  divergent  integral — if  an  integral  diverges,  then  we  should  conclude  that  we’ve  done 
something  wrong,  and  must  go  back  to  study  the  analysis  in  more  detail. 

In  contrast,  the  approach  taken  in  this  appendix  never  invokes  a  principal  value  to 
“fix”  a  divergent  integral,  because  it  has  no  divergent  integrals.  Rather,  the  principal 
value  was  deliberately  introduced  in  (A20)  as  a  way  of  joining  cos  cur  to  sin  cur  to  create 
a  complex  exponential  that  would  bring  the  full  power  of  Fourier  theory  to  bear,  enabling 
us  to  convert  sgnu;  to  a  Fourier  integral  in  (A27);  this  then  allowed  the  one-sided  Fourier 
integral  (A7)  to  be  calculated.  There  is  nothing  unusual  about  this  use  of  a  generalised 
function  to  perform  a  difficult  integral;  anyone  who  is  familiar  with  the  delta  function  5(t) 
in  Fourier  theory  should  have  no  trouble  accepting  that  it  has  a  sister  function,  cr(f).  In 
fact,  setting  S(co)  =  1  in  (A35)  gives 

/  e1^  dw  =  -  [1  +  i#l  2vr 6(t) 

Jo  2 

=  7T 5(t)  +  iu(t)/ir  *  nd(t) 

=  nd(t)  +  icr(t) ,  (A36) 

showing  that  the  two  generalised  functions  appear  side  by  side  in  this  fundamental  one¬ 
sided  Fourier  integral.12  Equation  (A36)  also  shows  that  cr(t)  is  the  Hilbert  transform 
of  7r<5(f);  again,  a  close  relationship. 

How  do  we  invert  the  Hilbert  transform?  This  can  be  done  through  the  evaluation 
of  a(t)  *  a(t ).  First,  write  the  Fourier  transform  (J  of  some  f(t)  as 

/OO 

e^fit)  d t  =  g(u) , 

-OO 

l  /*°° 

so  that  =  —  /  e~xwtg(u)du  =  f(t ) .  (A37) 

Now,  using  the  idea  that  the  Fourier  transform  of  a  convolution  equals  the  product  of  the 
individual  Fourier  transforms,  consider 

JF{a(t)  *  a(t)}  =  JF{a{t)}  x  JF {a{t)}  =  (ivrsgnu;)2  =  {  ^  ^  ^  °  (A38) 

[  0  oj  =  0  . 

Next,  invoke  the  inverse  Fourier  transform  of  both  sides  of  (A38),  noting  that  the  remov¬ 
able  discontinuity  at  u  =  0  can  be  ignored  for  the  purpose  of  doing  the  integral: 

a(t)  *  a(t)  =  9r~1{— 7T2}  ==  — —  [  e~lujt  duj  =  —TT25(t).  (A39) 

2-7T  7-oo 

(This  shows  that  a(t)/ir  is  a  “convolutional  square  root  of  —1”,  in  the  sense  that  5(t )  is 
to  convolution  what  1  is  to  multiplication:  convolving  any  function  with  5(t )  leaves  that 
function  unchanged.)  Now,  the  identity  cr(t)  *  a (t)  =  —  TT2S(t)  enables  us  to  invert  the 
Hilbert  transform  easily.  Consider  the  following,  with  “(i)”  omitted  for  clarity: 

s  =  9F{s}  =  s  *  a /it ,  (A40) 

12The  identity  (A36)  was  established  by  an  entirely  different  route  that  used  no  Fourier  transforms  in 
Chapter  11  of  [4]. 
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so  that 


9fi{s}  =  s*cr/ir  =  s  *  a /n  *  cr /n 

=  s  *  —8  =  — s  =  .  (A41) 

We  conclude  that  !H~l  =  —9fi.  Note  also  from  this  analysis  that 

9-[ {8}  =  cr / tt  ,  Jfi{a/Tr}  =  —8,  (A42) 

which  shows  more  of  how  the  Hilbert  transform  connects  these  two  generalised  functions. 


A. 3  Numerical  Evaluation  of  the  Hilbert  Transform 


We  can  express  the  Hilbert  transform  as  an  integral  over  time  in  parallel  streams  here  that 
swap  the  roles  of  s  and  a,  because  convolution  is  commutative: 


s(t)  =  cr(f)/7r  *  s(t) 


1 

7 r 

1 

7 r 


-  r)  dr 
dr 


=  s(t)  *  cr(f ) / 7T 

—  /  s(r)  a(t  —  r)  dr 

7T  J — oo 

s(r) 


=  1  i 


—  OO 
OO 


7T 


J—  OO  t 


dr . 


(A43) 


These  integrals  are  useful  in  theoretical  work  to  calculate  the  analytic  signal  s  +  is  for 
a  simple  real  signal  s(t).  But  they  can  be  difficult  to  apply  numerically  when  s(t)  is 
more  complicated.  When  transforming  a  signal  from  samples  that  have  all  been  collected 
(i.e.  we  are  not  examining  the  problem  of  how  best  to  transform  in  real  time  here),  an  easier 
approach  is  to  Fourier-transform  any  of  these  integrals  in  (A43),  recognising  that  they  are 
convolutions,  which  converts  them  to  a  product  of  Fourier  integrals,  which  is  then  inverse 
Fourier  transformed.  This  procedure  can  use  the  discrete  Fourier  transform  and  its  inverse, 
and  is,  in  fact,  completely  identical  to  implementing  the  frequency-domain  procedure  that 
we  used  to  define  the  analytic  signal  in  the  paragraph  immediately  preceding  (A6)  above. 

To  reiterate,  to  form  the  analytic  signal  sc(t)  of  a  real  signal  s(t),  first  Fourier- 
transform  s(t)  to  produce  its  spectrum.  Then  set  the  weights  of  all  negative  frequencies  to 
zero,  double  the  weights  of  all  positive  frequencies,  and,  to  provide  a  sense  of  Fourier  con¬ 
tinuity,  leave  the  weight  of  the  zero  frequency  unchanged:  this  is  equivalent  to  multiplying 
the  spectrum  by  0(uj)  using  (A15).  Finally  inverse  Fourier-transform  this  new  spectrum 
to  produce  sc(t). 

For  a  signal  composed  of  a  modulated  sinusoid,  the  value  of  the  modulation  (i.e.  en¬ 
velope)  at  any  time  can  be  defined  as  the  length  of  the  phasor  representing  the  complex 
number  sc(t): 

envelope  at  time  t  =  |sc(t)|  =  \J s2{t)  +  s2(t)  .  (A44) 

The  phase  <f>(t)  can  be  defined  as  the  argument  of  the  complex  number  sc(t): 

cos  fi(t)  =  ,  sin  cj)(t )  =  .  (A45) 

Mbl  Mb  I 

As  referred  to  in  (7.6),  at  any  moment,  the  instantaneous  or  “currently  dominant”  angular 
frequency  is  approximately  (j>'(t). 
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Sampling  at  200  Hz 


Figure  Al:  An  example  of  using  the  analytic  signal  to  construct  the  enve¬ 
lope  (black)  of  the  real  signal  in  (A46)  (red),  together  with  the  dominant  or 
running  frequency  at  each  time  step  (green,  using  the  same  y- axis  for  conve¬ 
nience  but  actually  different  units).  The  blue  curve  is  the  Hilbert  transform 
of  the  red  curve:  we  get  a  sense  of  the  90°  lag  referred  to  just  after  (A17). 
The  200  Hz  sampling  rate  gives  a  Nyquist  frequency  of  1 00  Hz,  more  than 
adequate  to  reveal  the  signal’s  structure. 


Examples  of  constructing  the  envelope  and  instantaneous  frequency  of  a  modulated 
signal  are  shown  in  Figures  Al  and  A2.  We  have  used  a  combination  of  two  modulation 
frequencies,  1  Hz  and  2  Hz,  on  a  carrier  with  frequency  20  Hz: 


s(t) 


(7  +  5  sin  27rf  sin  2n2t)  x  ^in  27r20ti  . 

'  .UYolop,'  '  Carrier 


(A46) 


For  Figure  Al  we  sample  at  200  Hz,  giving  a  100  Hz  Nyquist  frequency  that  is  of  course 
well  able  to  detect  the  20  Hz  carrier.  The  red  curve  in  the  plot  is  s(t )  and  the  blue  curve 
is  s(t),  which  is  seen  to  lag  the  red  curve  by  about  90°,  consistent  with  the  comment  just 
after  (A17).  The  black  envelope  is  calculated  as  the  modulus  of  sc(t)  (A44)  at  each  time, 
and  clearly  fits  the  signal  (red)  quite  well.  The  phase  is  calculated  from  (A45)  and  then 
numerically  differentiated:  the  phase  increment  at  each  time  step  is  divided  by  the  step’s 
time  increment  to  give  the  dominant  frequency  at  that  time,  which  is  plotted  in  green 
(using  the  same  y-axis  for  convenience  but  actually  different  units).  As  can  be  seen,  this 
dominant  frequency  matches  that  of  the  20  Hz  carrier  very  closely. 

Figure  A2  shows  results  for  various  sampling  rates  that  progressively  reduce  to  well 
below  carrier  Nyquist.  Although  the  carrier  is  generally  not  well  detected,  the  modulation 
(black  curve)  is  certainly  well  estimated,  showing  that  the  Hilbert  transform  is  quite  robust 
to  lower-than-Nyquist  sampling  rates. 
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Sampling  at  42  Hz 


time  (s) 

Sampling  at  20  Hz 


time  (s) 

Sampling  at  1 1  Hz 


time  (s) 


Sampling  at  37  Hz 


time  (s) 

Sampling  at  19  Hz 


time  (s) 

Sampling  at  7  Hz 


time  (s) 


Figure  A2:  Results  for  the  signal  of  Figure  Al,  but  now  sampled  at  progres¬ 
sively  lower  frequencies.  The  labels  and  legends  of  these  plots  are  the  same  as 
those  of  Figure  A 1  but  are  mostly  omitted  for  clarity.  Sampling  at  multiples 
of  the  carrier  (such  as  the  20  Hz  at  mid  left)  gives  no  useful  result. 
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Appendix  B  Tutorial  on  Correlation  and 

Convolution 

The  process  of  correlating  two  sequences  of  numbers  is  central  to  much  of  this  report,  and 
a  good  understanding  of  it  is  necessary  to  follow  the  discussions  that  accompany  equations 
such  as  (4.6).  This  appendix  gives  step-by-step  descriptions  of  the  processes  of  correlation 
and  convolution.  The  two  procedures  are  almost  identical,  so  that  a  knowledge  of  one  is 
easily  transferred  to  the  other.  Correlation  (“*”)  detects  the  extent  to  which  a  sequence 
of  numbers  (the  signal)  is  present  inside  a  possibly  longer  sequence  that  is  suspected 
to  contain  the  signal  together  with  noise.  In  contrast,  convolution  (“*”)  is  a  procedure 
that  calculates  the  output  of  a  linear  time-shift-invariant  system  given  some  input,  and 
is  heavily  involved  in  Fourier  transform  theory.  Rewriting  a  correlation  as  a  convolution 
enables  us  to  search  for  one  signal  inside  another,  by  using  a  standard  technique  for 
convolving  quickly  and  efficiently  that  is  based  on  the  discrete  Fourier  transform. 

Although  correlation  is  defined  as  essentially  convolution  below,  there  is  no  reason 
that  it  has  to  be;  convolution  is  simply  the  standard  “matched  filter”  choice  of  correlation 
procedure,  chosen  because  convolution  turns  out  to  maximise  the  ratio  of  what  results 
from  correlating  something  with  a  signal,  compared  to  what  results  from  correlating  with 
the  noise  that  is  invariably  attached  to  that  signal. 

B.l  The  Correlation  Procedure 

Consider  searching  for  a  sequence  [— 1,  2,  3],  which  constitutes  a  signal,  within  a  sequence 
that  might  contain  the  signal  together  with  some  noise.  We’ll  take  this  sequence  to  be 
[0.1,  0.2,  -1,  2.1,  3,  0.1]. 

Table  B1  shows  the  steps  of  the  procedure.  It  consists  of  moving  [—1,  2,  3]  in  from  the 
left  to  overlap  [0.1,  0.2,  —1,  2.1,  3,  0.1],  one  element  at  a  time,  and  at  each  step  summing 
the  products  of  pairs  of  numbers  that  overlap:  these  pairs  are  shown  in  red.  At  each  step, 
this  sum  is  output  as  the  next  element  of  the  correlation  sequence.  The  process  is  identical 
to  finding  the  euclidean  dot  product,  at  each  step,  of  the  two  currently  overlapping  arrays 
of  red  numbers.  We  write  the  correlation  as 

[-1,  2,  3]  *  [0.1,  0.2,  -1,  2.1,  3,  0.1]  = 

[0.3,  0.8,  -2.7,  4.1,  14.2,  4.2,  -2.8,  -0.1]  .  (Bl) 

The  numbers  to  be  correlated  will  generally  be  complex,  and  a  complex  conjugate 
must  be  used  in  the  correlation  process.  The  reason  was  given  on  page  17:  only  by  using 
a  complex  conjugate  will  this  procedure  give  a  meaningful  result:  a  peak  when  the  signal 
matches  up  with  itself.  Note  that  by  convention,  no  explicit  conjugation  is  written  with 
the  notation.  For  example  when  correlating  [5,  6+21,  3+8i]  with  [1+i,  3,  5— 4i,  7], 
we  write  the  result  as 

[5,  6+2i,  3+8i]  *  [1+i,  3,  5-4i,  7]  ,  (B2) 

but  the  correlation  is  done  by  moving  the  conjugated  sequence  [5,  6  — 2i,  3  —  81]  in  from 
the  left  over  [1+i,  3,  5  — 4i,  7]. 

The  correlation  in  (Bl)  has  a  maximum  14.2  as  its  fifth  element.  As  can  be  seen  from 
the  steps  of  the  procedure,  it  follows  that  the  best  estimate  of  where  the  signal  sequence 
begins  within  the  signal-plus-noise  sequence  is  the  third  entry  of  the  latter.  Locating  the 
best  estimate  like  this  is  the  idea  behind  the  analysis  of  Section  4.2. 
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Table  Bl:  The  array  [— 1,  2,  3]  is  correlated  with  [0.1,  0.2,  —1,  2.1,  3,  0.1] 
through  the  following  steps.  The  result  is  the  bottom  row  of  numbers 
[0.3,  0.8,  —2.7,  4.1,  14.2,  4.2,  —2.8,  —0.1].  At  each  step,  the  number  pairs 
to  be  multiplied  are  written  in  red. 

Correlation  step  1:  —1  2 

3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

Correlation  step  2:  —1 

2 

3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

0.8 

Correlation  step  3: 

-1 

2 

3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

0.8 

-2.7 

Correlation  step  4: 

-1 

2 

3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

0.8 

-2.7 

4.1 

Correlation  step  5: 

-1 

2 

3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

0.8 

-2.7 

4.1 

14.2 

Correlation  step  6: 

-1 

2 

3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

0.8 

-2.7 

4.1 

14.2 

4.2 

Correlation  step  7: 

-1 

2  3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

0.3 

0.8 

-2.7 

4.1 

14.2 

4.2  -2.8 

Correlation  step  8: 

-12  3 

0.1 

0.2 

-1 

2.1 

3 

0.1 

The  final  result: 

0.3 

0.8 

-2.7 

4.1 

14.2 

4.2  -2.8  -0.1 
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B.2  The  Convolution  Procedure 

The  series  of  multiplications  that  implement  a  correlation  (“*”)  is  in  fact  just  a  re-ordered 
version  of  the  set  that  implement  a  convolution  (“*”).  Convolution  is  simply  polynomial 
multiplication,  and  the  way  it  relates  to  correlation  is  via  (4.3): 

A*B  =  *  B.  (B3) 

There  is  no  mystery  here;  correlation  and  convolution  are  just  two  different  procedures  that 
happen  to  be  related  by  reversing  and  conjugating  one  of  the  sequences.  Equation  (B3) 
expresses  the  correlation  in  (Bl)  with  a  convolution: 

[-1,  2,  3]  *  [0.1,  0.2,  -1,  2.1,  3,  0.1]  =  [3,  2,  -1]  *  [0.1,  0.2,  -1,  2.1,  3,  0.1].  (B4) 

Table  B2  shows  the  steps  followed  in  the  convolution  on  the  right-hand  side  of  (B4).  The 
result  is 


[3,  2,  -1]  *  [0.1,  0.2,  -1,  2.1,  3,  0.1]  = 

[0.3,  0.8,  -2.7,  4.1,  14.2,  4.2,  -2.8,  -0.1]  ,  (B5) 

just  as  was  obtained  for  the  correlation  in  (Bl). 

We  are  working  with  real  numbers  only  for  convenience,  but  as  an  example  of  using  complex 
numbers,  (B2)  becomes 

[5,  6+2i,  3+8i]  *  [1+i,  3,  5-4i,  7]  =  [3-81,  6-21,  5]  *  [1+1,  3,  5-41,  7] 

=  [11-51,  17-201,  6-531,  58-901,  67-341,  35].  (B6) 

The  numbers  in  (B5)  are,  by  design,  the  coefficients  in  the  following  polynomial  mul¬ 
tiplication: 

(3  +  2 2  -  z2)  (0.1  +  0.2 z  -  z2  +  2 .lx3  +  3z4  +  O.lz5)  = 

0.3  +  0.8 z  -  2.7 z2  +  4.1  z3  +  14.2 z4  +  4.2 z5  -  2.8z6  -  0.1  z7  .  (B7) 

These  polynomials  are  defined  to  be  the  z-transforms  of  the  original  sequences  (z  is  re¬ 
placed  by  1/z  in  an  alternative  definition  of  the  transform).  This  means  that  convolving 
two  sequences  is  equivalent  to  multiplying  their  z-transforms  and  inverse  z-transforming 
the  result,  which  converts  the  polynomial  product  back  to  a  sequence  of  numbers — its 
coefficients.  This  procedure  forms  the  Convolution  Theorem  for  the  z-transform.  The 
polynomial  formed  from  A  =  {A0,  A1:  A2,  ■  ■  ■  }  is  the  z-transform  of  A: 

Z>{A}  =  A0  +  A1z  +  A2z2  +  ...  (B8) 

and  similarly  for  B.  The  procedure  in  (B7)  becomes 

Z{A }  x  Z{B}  =  Z{A  *  B}  ,  (B9) 

which  is  the  Convolution  Theorem.  This  correspondence  between  convolution  and  multi¬ 
plication  under  some  transform  is  a  central  theme  of  convolution  theory. 

The  z  that  appears  in  the  z-transform’s  polynomial  is  not  normally  given  any  value; 
it  serves  solely  to  convert  a  sequence  to  a  function.  But  if  we  set  z  equal  to  10  in  the 
polynomials  of  (B7),  the  convolution  becomes  a  straightforward  base-10  multiplication — 
except  that  the  digits  in  a  base-10  multiplication  are  usually  positive. 
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Table  B2:  The  array  [3,  2,  —1]  is  convolved  with  [0.1,  0.2,  —1,  2.1,  3,  0.1] 
by  the  procedure  below.  The  result  of  each  step  is  shown  in  red  under 
the  solid  line  for  that  step,  so  that  the  last  line  becomes  the  convolution: 

[0.3,  0.8,  —2.7,  4.1,  14.2,  4.2,  —2.8,  —0.1].  All  numbers  that  are  multiplied 
above  the  line — in  pairs,  from  outside  inwards — are  written  in  red.  We  ap¬ 
pend  2  zeroes  to  the  right-hand  sequence  ( 2  being  one  less  than  the  length  of 
the  left-hand  sequence);  the  only  function  of  these  is  to  facilitate  the  pairwise- 
multiplication  procedure  when  using  pen  and  paper. 

Convolution  step  1:  3  2  —1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

(= 

3  x  0.1) 

Convolution  step  2:  3  2  —1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

0.8 

(=3x0.2 +  2x0.1) 

Convolution  step  3:  3  2  —  1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

0.8 

-2.7 

(= 

3  x  —1  +  2  x  0.2  - 

1  x  0.1) 

Convolution  step  4:  3  2  —  1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

0.8 

-2.7 

4.1 

(=  3  x  2.1  +  etc. 

) 

Convolution  step  5:  3  2  —  1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

0.8 

-2.7 

4.1 

14.2 

Convolution  step  6:  3  2  —  1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

0.8 

-2.7 

4.1 

14.2  4.2 

Convolution  step  7:  3  2  —1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

0.3 

0.8 

-2.7 

4.1 

14.2  4.2 

-2.8 

Convolution  step  8:  3  2  —  1  * 

0.1 

0.2 

-1 

2.1 

3  0.1 

0 

0 

The  final  result: 

0.3 

0.8 

-2.7 

4.1 

14.2  4.2 

-2.8 

-0.1 
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In  fact,  the  digits  needn’t  be  positive;  there  is  nothing  wrong  with  writing  e.g.  28  as  32,  meaning 
30  —  2,  so  that  a  multiplication  such  as  54  x  28  becomes  54  x  32  =  1528  =  1512.  Multiplication 
is  often  done  this  way  by  fast  human  mental  calculators,  because  it  converts  the  “difficult”  mul¬ 
tiplication  by  8  into  an  easy  multiplication  by  2,  with  the  only  price  being  a  minor  amount  of 
subtraction  needed. 

The  carrying  step  necessary  for  base-10  multiplication  is  not  a  part  of  the  convolution  of 
course;  rather,  it’s  simply  a  consequence  of  the  fact  that  any  base-n  system  uses  single  digits  to 
represent  numbers  less  than  n. 

Conventionally  of  course,  the  digits  making  up  a  base-10  number  are  written  in  the  reverse 
direction  to  that  of  the  coefficients  in  (B7),  so  that  the  convolution  that  comprises  a  base-10 
multiplication  ends  up  being  done  from  right  to  left.  You  can  see  this  in  more  detail  when 
multiplying  1234  by  567,  following  the  steps  in  Table  B3.  We  first  express  the  numbers  as 
sequences  in  base-10  by  z-transforming  them  (with,  effectively,  z  set  equal  to  10,  although 
we  are  so  familiar  with  mentally  equating  a  number  with  its  “10-transform”  that  we  might 
not  be  explicitly  aware  of  doing  this  step).  Then  we  convolve  the  sequences,  and  finally 
we  inverse  z-transform  the  resulting  sequence  back  to  a  number — again,  a  step  that  we 
might  not  explicitly  be  aware  of  doing,  since  we  take  for  granted  the  actually  deep  notion 
of  equality  between  a  number  and  its  base-10  representation.  In  other  words,  we  multiply 
numbers — in  principle  a  difficult  task  that  involves  lots  of  counting — by  convolving  their 
“10-transforms”,  which  is  a  much  easier  task.  We  see  then,  that  what  we’re  accustomed 
to  think  of  as  multiplying  two  numbers  written  with  digits  is  in  fact  a  convolution  of  those 
digits — with  added  “carrying  steps”  that  are  necessary  when  using  a  finite  set  of  digits  to 
represent  the  numbers. 

Convolution  and  the  Discrete  Fourier  Transform 

Equation  (B9)  holds  fairly  trivially.  A  similar  expression  results  when  the  z-transform  is 
replaced  by  the  discrete  Fourier  transform  2),  which  transforms  a  sequence  of  complex 
numbers  to  another  sequence  of  the  same  length.  If  the  sequences  A  and  B  have  the  same 
length,  then  the  Convolution  Theorem  states 

©{A}  x  D{B}  =  2 1{A  *  B}  ,  (BIO) 

where  the  multiplication  “  x  ”  of  sequences  P  and  Q  is  defined  pointwise:  ( P  x  Q)n  =  PnQn. 
When  A  and  B  have  different  lengths,  a  slight  alteration  to  (BIO)  is  needed.  Their  con¬ 
volution  A*  B  has  length  L  =  length(A)  +  length (i?)  —  1.  Append  zeroes  to  A  and  B 
to  make  new  arrays  A'  and  B'  respectively  that  both  have  length  L.  (Note,  this  is  not 
related  to  the  mostly  outdated  practice  of  zero  padding,  discussed  ahead  in  Appendix  D.) 
With  the  definition  of  the  DFT  used  in  this  report  [equations  (6.1),  (6.2),  and  Table  Cl 
in  Appendix  C],  the  Convolution  Theorem  becomes 

LtD{A'}  x  £>{B'}  =  P{A  *  B}  .  (Bll) 

The  Matlab  functions  f  f  t  and  if  ft  differ  in  their  normalisation  from  my  Df  t  and  Df  tlnverse 
of  Table  Cl;  Matlab ’s  normalisation  makes  (Bll)  slightly  simpler: 

®Matlab{“4/}  X  ‘Dua.tlabiB' }  =  ®Matlab{^  *  B}  .  (B12) 

On  the  other  hand,  f  f  t  doesn’t  immediately  give  the  correct  peak  heights  to  spectra  made 
from  it,  whereas  Dft  does.  Always  check  the  normalisation  used  in  your  own  choice  of 
Fast  Fourier  Transform  routine. 
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Table  B3:  Multiplying  1234  x  567  to  get  699,678  by  recognising  that  the 
process  is  a  convolution  plus  the  mechanism  of  “carrying”  which  is  necessary 
to  represent  the  result  in  a  base- 10  system.  We  prepend  the  left-hand  factor 
with  a  number  of  zeroes  that  is  one  less  than  the  number  of  digits  in  the 
right-hand  factor.  This  procedure  is  the  mirror  image  of  the  procedure  in 
Table  B2. 


Convolution  step  1:  0 

0 

1 

2 

3 

4 

X 

5 

6 

7 

4x7  = 

28: 

28 

Convolution  step  2:  0 

0 

1 

2 

3 

4 

X 

5 

6 

7 

3x7+4x6+  carried  2  = 

47: 

47 

8 

Convolution  step  3:  0 

0 

1 

2 

3 

4 

X 

5 

6 

7 

2x7  +  3x6  +  4x5  +  carried  4  = 

56: 

56 

7 

8 

Convolution  step  4:  0 

0 

1 

2 

3 

4 

X 

5 

6 

7 

1x7  +  etc.  = 

39: 

39 

6 

7 

8 

Convolution  step  5:  0 

0 

1 

2 

3 

4 

X 

5 

6 

7 

0x7  +  etc.  =  19: 

*9 

9 

6 

7 

8 

Convolution  step  6:  0 

0 

1 

2 

3 

4 

X 

5 

6 

7 

The  final  answer:  6 

9 

9 

6 

7 

8 

Equation  (Bll)  is  a  highly  efficient  method  for  convolving  A  and  B.  For  example,  we 
can  use  it  to  reproduce  the  result  from  Table  B2  with  this  Matlab  code: 

A  =  [3  2  -1] ; 

B  =  [0.1  0.2  -1  2.1  3  0.1]; 

A_prime  =  [A  zeros (1 , length (B) -1) ] ; 

B_prime  =  [B  zeros ( 1 . length ( A )- 1 ) 1 : 

( length (A)+length (B) -1)  *  DftlnverseC  Dft(A_prime)  .*  Dft (B_prime 

7,  The  above  line  is  identical  to 
if f t (  f f t ( A  prime)  .*  fft(B  prime)  ) 

We  get  the  expected  result  of 

[0.3  0.8  -2.7  4.1  14.2  4.2  -2.8  -0.1] 

For  very  long  arrays,  convolution  never  uses  the  expensive  procedure  of  Table  B2.  Instead, 
it  always  uses  the  Convolution  Theorem,  specifically  (Bll)  or  (B12).  But  understanding 
the  procedure  of  Table  B2  will  make  you  a  very  adept  convolver. 

On  a  final  note,  this  replacing  of  a  difficult  convolution  by  an  easy  pointwise  multipli¬ 
cation  via  the  discrete  Fourier  transform  is  the  reverse  idea  to  the  multiplication  example 
in  Table  B3.  When  learning  arithmetic,  we  effectively  replace  difficult  multiplication  by 
easy  convolution  via  a  “10-transform” — although  few  would  ever  describe  it  that  way. 
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Appendix  C  Sample  Discrete  Fourier  Transform 

Code 

Here  we  list  two  examples  of  Matlab  and  Mathematica  code  to  implement  a  DFT.  Our 
convention  is  to  plot  only  frequencies  with  absolute  value  less  than  the  Nyquist  frequency; 
thus  negative  frequencies  are  plotted  as  negative  frequencies,  and  not  wrapped  around 
to  become  positive  frequencies  greater  than  the  Nyquist  frequency.  In  particular,  our 
convention  allows  a  plot  of  the  DFT  of  a  gaussian  function  to  be  another  gaussian.  This 
mimics  the  fact  that  the  continuous  Fourier  transform  of  a  gaussian  function  is  another 
gaussian. 

Table  Cl  shows  an  example  of  Matlab  code  that  implements  the  DFT  in  (6.1)  and 
inverse  DFT  in  (6.2).  It  might  appear  that  the  code  here  for  a  matrix  could  also  be  used 
for  an  array,  thus  eliminating  the  need  for  the  “if-else”  block.  But  that  isn’t  so,  because 
the  “,  1”  in  the  fftshift  command  is  really  only  appropriate  for  matrices;  when  used  on 
an  array,  it  gives  a  different  result  for  a  row  array  than  for  a  column  array. 

Table  C2  has  an  example  of  Mathematica  code  that  does  the  same  job. 


Table  Cl:  An  example  of  Matlab  code  that  implements  the  DFT  in  (6.1) 
and  inverse  DFT  in  (6.2).  See  the  text  above  for  the  use  of  1”  in  the 
fftshift  command. 

The  DFT  can  be  implemented  in  Matlab  with 
function  X  =  Dft(x) 

7,  If  x  is  an  array: 
if  isvector (x)  ==  1 

X  =  fftshift(fft(x))/length(x)  : 
else 

'/,  x  is  a  matrix,  so  act  on  each  column  separately. 

X  =  fftshift(fft(x)  ,l)/size(x,l); 
end 


The  inverse  DFT  can  be  implemented  with 
function  x  =  Df tlnverse (X) 

%  If  X  is  an  array: 
if  isvector (X)  ==  1 

x  =  ifft(ifftshift(X))  *  length (X) : 
else 

7,  X  is  a  matrix,  so  act  on  each  column  separately, 
x  =  ifft(ifftshift(X,l))  *  size (X  ,  1)  ; 
end 
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Table  C2:  An  example  of  Mathematica  code  that  implements  the  DFT 
in  (6.1)  and  inverse  DFT  in  (6.2) 

(*  These  do  what  Matlab's  fftshift  and  ifftshift  do. 

For  x  any  list  , 

Ff t shif t [ I f f t shi f t [x] ]  =  If ftshift [Fftshift  [x] ]  =  x.  *) 

Ff tshif t [x_List]  : =  (  NN  =  Length  [xl  : 

FlattenTlTake  Tx.  -Floor  f  NN  /  2 1  1  ,  Take  T  x  .  _Cei_l_ing  [NN  /  2  ]  ]  }]  ) 

If ftshift [x_List]  : =  (  NN  =  Length [x] ; 

Flatten  [{Take  [x  ,  -Ceiling  [NN/2] ]  ,  Take [x ,  Floor  [NN / 2 ] ] }]  ) 

Begin ["Private'"] 

(*  These  do  what  my  Matlab  functions  Dft  and  Dftlnverse  do, 
but  just  on  an  array.  They  are  defined  them  here  solely 
for  use  in  building  Dft  and  Dftlnverse  in  a  moment. 

For  x  any  list,  Df tlnverseVector [Df tVector [x] ]  = 

Df tVector [Df tlnverseVector  [x] ]  =  x.*) 

Df tVector [x_List ]  :=  Fftshift T Fourier  [x. 

FourierParameters  ->  {-1,  -1}]] 

Df t Inver s eVect or [X_List ]  :=  InverseFourier [Ifftshift [X] , 

FourierParameters  ->  {-1,  -1}] 

End  n  (*  "Private'"  *) 


(*  These  do  what  my  Matlab  functions  Dft  and  Dftlnverse  do, 
on  arrays  and  on  matrices.  *) 

Dft [x_List]  :  = 

If [Length [Dimensions  Txll  ==  1, 

(*  x  is  an  array  *) 

Private ' Df tVector  [x]  , 

(*  else  x  is  a  matrix  *) 

Transpose [ 

Table  T Private ' DftVector  [  x  [  [  All  .ill  1  ,  {i,  Length  [  x[[l]]  ] }] 

1 

1 

Df t Inver se [X_Li st ]  := 

If  T Length  T Dimensions [Xl 1  ==  1, 

(*  X  is  an  array  *) 

Private ' Df tlnverseVector  [X]  , 

(*  else  X  is  a  matrix  *) 

Transpose [ 

Table  F Private ' DftlnverseVector  [  X  [  [ All  .ill  ], 

{i,  Length  [  X  [  [  1 ] ]  ] >] 

1 

1 
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Appendix  D  Zero  Padding  in  the  Discrete 

Fourier  Transform 

The  appending  with  zeroes  that  was  required  in  the  analysis  around  (Bll)  to  convolve 
two  sequences  with  the  DFT  is  an  exact  procedure:  it  gives  the  exact  same  result  that’s 
obtained  when  we  convolve  the  two  sequences  using  the  approach  of  Table  B2.  This 
appending  is  distinct  from  the  historical  zero  padding  of  data  mentioned  in  Section  6. 
Older  implementations  of  the  FFT  required  the  length  of  the  data  set  to  be  a  power  of  2, 
and  so  a  data  set  with  length  other  than  this  would  have  zeroes  appended  to  create  one 
whose  length  was  a  power  of  2.  This  zero  padding  changed  the  data  and  so  introduced 
spurious  frequencies  to  the  spectrum.  These  unwanted  frequencies  were  tolerated  because 
only  by  zero  padding  to  a  power-of-2  length  could  the  FFT  be  used  at  all.  But  current 
FFT  algorithms  are  very  fast  even  when  the  data  length  is  not  a  power  of  2,  so  this  zero 
padding  is  no  longer  necessary. 

To  see  why  zero  padding  changes  a  Fourier  transform,  recall  the  central  fact  that  the 
DFT  is  a  set  of  amplitudes  of  sinusoids.  Because  sinusoids  are  periodic,  the  DFT  effectively 
“assumes”  that  the  sequence  of  numbers  being  transformed  is  periodic.  The  result  is  that 
the  corresponding  Fourier  series  replicates  the  original  sequence  endlessly  forward  and 
backward  in  time.  For  example,  the  5-element  sequence  [1,2, 3,4,  5]  is  interpreted  by  the 
DFT  to  be  the  infinite  data  set 

...1,2, 3, 4, 5,1, 2, 3, 4, 5,1, 2, 3, 4, 5,...  (Dl) 

What  we  mean  here  is  that  the  Fourier  series  that  returns  the  numbers  [1,2, 3, 4,  5]  as 
a  function  of  times  equal  to,  say,  1,2, ...  ,5  will  be  the  same  as  the  series  that  returns 
the  numbers  [1,  2,  3, 4,  5, 1,  2, 3, 4,  5]  as  a  function  of  times  1,2,...,  10.  The  DFTs  of  these 
two  sequences  differ  only  by  the  presence  of  interspersed  zeroes  in  the  DFT  of  the  second 
sequence;  but  those  zeroes  are  weights  of  sinusoids,  and  so  the  end  result  is  that  both  se¬ 
quences  give  rise  to  the  same  non-zero  sinusoids:  the  same  Fourier  series.  (Verify  this  using 
Dft — but  not  fft,  which  doesn’t  normalise  appropriately.)  But  zero  padding  [1,2,3, 4,  5] 
to  the  next  power-of-2  length  alters  it  to  [1,  2, 3, 4,  5,  0, 0,  0] ,  which  is  then  interpreted  by 
the  DFT  as 


...  1,  2,  3, 4,  5, 0,  0, 0, 1,  2,  3, 4,  5,  0, 0, 0, 1,  2, 3, 4,  5,  0,0,0,... 


(D2) 


The  spectrum  of  this  sequence  is  certainly  different  to  that  of  (Dl).  In  essence,  the 
triplets  of  zeroes  in  (D2)  are  treated  as  a  constant  stream  of  data  that  requires  compar¬ 
atively  many  frequencies  to  generate  it,  because  many  sinusoids  are  required  to  cancel 
each  other  to  produce  a  “zero-slope”  sequence.  The  spectra  generated  from  [1,2,3, 4,  5] 
and  [1,  2, 3, 4,  5,  0, 0,  0]  are  shown  in  Figure  Dl  on  the  following  page.  See  the  end  of  this 
appendix  for  some  discussion  of  the  right-hand  plot. 

Some  users  of  the  DFT  object  that  the  DFT  algorithm  effectively  replicates  their  data 
into  the  past  and  future.  They  wish  to  let  values  of  zero  represent  unknown  values  of  the 
system  under  observation,  and  so  will  append  many  zeroes  to  their  data  before  Fourier 
transforming  it.  The  length  of  the  padded  sequence  need  not  be  a  power  of  2,  so  this  zero 
padding  is  not  related  to  the  historical  zero  padding  required  by  early  FFT  algorithms. 
When  this  zero-padded  sequence  is  effectively  replicated  by  the  DFT,  the  zeroes  have  the 
effect  of  pushing  the  data  replications  very  far  into  the  past  and  future  from  the  sequence 
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Frequency  (Hz) 


Figure  Dl:  Left:  Spectrum  of  [1,2, 3, 4, 5],  sampled  once  per  second. 
Right:  Spectrum  of  [1, 2,  3, 4,  5, 0, 0,  0] ,  sampled  once  per  second.  More  and 
different  frequencies  are  necessary  to  construct  the  Fourier  series  that  repre¬ 
sents  the  second  data  set. 


that  is  “located  near  the  present”,  so  to  speak.  But  there  is  no  reason  why  this  zero 
estimate  of  what  was  never  measured  should  be  any  more  realistic  than  the  replicating 
alternative.  Indeed,  if  we  really  don’t  want  the  Fourier  transform  to  treat  our  data  as 
periodic,  then  perhaps  we  shouldn’t  be  using  a  Fourier  transform  in  the  first  place. 

Nevertheless,  suppose  we  append  1000  zeroes  to  [1,2, 3, 4,  5]  and  then  transform  the 
resulting  long  array.  The  resulting  spectrum  is  shown  in  Figure  D2.  (The  high  density 
of  plot  points  in  that  figure  obliges  us  to  omit  the  stems  that  were  used  to  highlight  the 
plot  points  of  Figure  Dl.)  The  amplitudes  being  plotted  are  now  quite  small.  This  is 
reasonable,  given  that  each  amplitude  is  the  amount  present  of  its  corresponding  sinusoid, 
and  the  more  sinusoids  we  have,  the  less  we’ll  generally  need  of  each. 

You  can  think  of  zero-padding  data  as  equivalent  to  multiplying  a  longer  set  of  data 
by  a  top-hat  function  that  zeroes  the  last  segment  of  the  data.  Fourier-transforming  this 
product  of  data  and  top  hat  is  identical  to  convolving  the  Fourier  transforms  of  the  data 
and  top  hat.  The  top  hat  transforms  to  a  sine  function,  so  the  effect  is  to  convolve  the 
correct  spectrum  with  a  sine.  This  smears  the  spectrum  out,  introducing  the  spurious 
frequencies  mentioned  at  the  start  of  this  appendix.  (But  remember  that  this  smearing 


Figure  D2:  The  DFT  of  the  array  [1,  2, 3, 4,  5,  0, ...  ,  0],  where  1000  zeroes 
are  appended  after  the  5.  The  spectrum  is  still  discrete,  but  the  “stems'’  that 
are  used  in  Figure  Dl  are  omitted  here. 
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has  no  relevance  if  we  are  using  the  FFT  as  a  fast  way  of  convolving  two  different-length 
arrays.) 

Zero  padding  is  sometimes  described  as  a  way  of  returning  a  finer  or  better  frequency 
resolution.  But  the  DFT  cannot  work  magic:  it  cannot  somehow  go  back  out  into  the  field 
and  return  a  whole  new  set  of  data  points.  Zero  padding  smoothens  the  spectrum  with 
a  sine  function  and  returns  a  denser  set  of  frequencies,  but  that  doesn’t  imply  that  these 
frequencies  actually  exist  in  the  data. 

You  can  find  zero  padding  described  in  the  literature  as  relatively  harmless  because 
the  zeroes  are  “not  new  data”  as  they  “don’t  contribute  to  any  sums  in  the  FFT”.  But, 
of  course,  the  FFT  doesn’t  know  which  of  its  input  numbers  are  data  and  which  are  not 
data;  it  simply  processes  an  array  that  it  can  only  treat  uniformly.  So  the  padded  zeroes 
effectively  become  new  data — but  false  data.  And  that  the  zeroes  don’t  contribute  to  any 
sums  is  quite  irrelevant.  The  FFT  requires  as  input  the  total  number  of  data  points,  and 
that  number  certainly  increases  when  zeroes  are  appended,  which  changes  the  FFT  output 
non-trivially.  Analogously,  if  we  were  to  form  a  sum  to  calculate  the  mean  of  that  data 
and  then  do  the  same  with  appended  zeroes  included,  the  zeroes  wouldn’t  contribute  to 
that  sum  either—  but  they  certainly  would  change  the  mean,  as  indeed  they  should. 

Another  way  of  seeing  that  zero  padding  has  no  a  priori  validity  in  building  a  spectrum 
is  to  consider  two  researchers  who  independently  collect  identical  data:  say  100  numbers. 
Their  FFT  algorithm  requires  its  data  set  to  have  length  a  power  of  2,  so  Researcher  A 
pads  with  28  zeroes.  Researcher  B,  on  the  other  hand,  collects  28  more  data  points  and 
is  astonished  to  find  that  they’re  all  zeroes.  A  and  B  will  end  up  Fourier  transforming 
identical  sets  of  128  numbers.  But  B  has  measured  something  that  A  has  not,  and  so  B 
should  surely  produce  a  spectrum  with  more  information  than  that  obtained  by  A,  and 
yet  does  not;  instead,  Researcher  A  has  by  sheer  coincidence  or  luck  produced  the  exact 
same  spectrum  obtained  by  B. 

The  acid  test  of  any  FFT  algorithm  is  “Sample  one  or  more  noiseless  sinusoids  in  the 
way  described  for  Figure  5  on  page  26,  and  then  plot  the  spectrum  produced  by  the  FFT. 
This  spectrum  should  be  non-zero  only  for  each  of  the  frequencies  used,  as  in  Figure  5”. 
The  reason  why  no  extra  frequencies  should  appear  is  because  that  procedure  samples  the 
sinusoids  in  such  a  way  that  when  the  DFT  effectively  strings  copies  of  the  data  together, 
what  results  still  represents  one  or  more  pure  sinusoids.  (More  generally  we  might  choose 
fractional  frequencies,  but  the  sampling  required  to  preserve  each  as  a  single  sinusoid  in 
the  replication  becomes  more  delicate  then.) 

For  example,  in  Matlab’s  notation,  choose  a  sampling  frequency  sf  that  is  a  natural 
number  (i.e.  anything  in  1,  2,  3, ... )  and  form  a  set  of  data  that  oscillates  at,  say,  4  Hz: 

data  =  sin  (2*pi  *  4  *  linspace (0 ,  1-1/sf,  sf ) ) 

where  the  sampling  frequency  is  above  Nyquist’s  8  Hz.  (The  sampling  frequency  must  be  a 
natural  number  in  this  case  to  ensure  that  Matlab  builds  an  array  of  samples  precisely  cor¬ 
responding  to  each  time  in  the  array  0  :  1/sf  :  1-1/sf.  That  is,  the  final  value  1-1/sf 
must  equal  a  natural  number  of  increments  1/sf.)  Yes,  the  sampling  required  by  this  test 
is  quite  artificial,  but  that  has  no  bearing  on  the  use  of  this  test  to  ascertain  whether  an 
FFT  algorithm  works  or  not.  The  FFT  of  data  should  take  on  the  value  of  exactly  0.5  at 
the  frequencies  of  exactly  4  and  —4,  and  be  exactly  zero  for  all  other  frequencies.  Zero¬ 
padding  data  before  FFT-ing  it  will  fail  this  test.  (Strangely,  some  zero-padders  will  insist 
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that  the  two  frequencies  here  of  exactly  4  and  —4  are  incorrect.  But  those  frequencies  are 
certainly  correct  because  they  represent  a  sine  function  of  frequency  exactly  4  Hz.) 

A  final  comment:  the  right-hand  spectrum  in  Figure  D1  includes  the  negative  Nyquist 
frequency  (—0.5)  but  not  the  positive  one.  Is  this  correct?  It  certainly  is;  the  DFT  has 
8  elements  here,  and  because  the  zero  frequency  is  always  present  for  both  even-  and 
odd-length  data  sets,  only  the  negative  Nyquist  frequency  will  be  represented  in  the  DFT 
when  the  data  has  even  length.  But  the  weighting  of  this  frequency  includes  the  weighting 
of  the  positive  Nyquist  frequency.  For  a  closer  inspection  of  this  figure,  use  the  DFT  (6.1) 
and  its  inverse  (6.2).  The  array  of  data  is  [1,  2,  3, 4,  5, 0,  0,  0] .  Call  these  data  elements  xk 
for  k  =  0  to  7  respectively.  The  DFT  of  this  array  is  [X_4,  X_3, . . . ,  X3],  whose  elements 
are  respectively 

0.375,  -0.3232  + 0.1553  i,  0.375  -0.251,  -0.6768  +  0.9053  i,  1.875, 

-0.6768  -0.90531,  0.375  +  0.251,  -0.3232  -  0.1553  i .  (D3) 

The  absolute  values  of  these  numbers  are  plotted  in  the  right-hand  spectrum  of  Figure  Dl, 
versus  their  frequencies,  which  are  the  set  oin/N  in  (6.1),  namely  —4/8,  —3/8,  —2/8, . . .  ,3/8. 
These  frequencies  now  insert  into  (6.2)  to  give  the  inverse  transform,  for  k  =  0  to  7: 

xk=  ]T  Xn  ei27rWJV 

n=— 4 

=  0.375  eiri27rfc  +  (-0.3232  +  0.1553i)  e^i2nk  +  (0.375  -  0.25  i)  e^ri2nk 
+  (-0.6768  +  0.9053 i)  e^i2nk  +  1.875  eli27rfc  +  (-0.6768  -  0.90531)  e^i2nk 
+  (0.375  +  0.25  i)  eii2?rfc  +  (-0.3232  -  0.1553 i)  e^i2nk.  (D4) 

Setting  k  =  0  to  7  here  returns  the  numbers  [x0,  Xi,  x2,  x3,  X4,  x5,  x6,  x7]  as  [1,  2, 3, 4,  5, 0,  0,  0] , 
as  expected. 
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Appendix  E  Calculating  a  Signal-to-Noise  Ratio 

Here  we  calculate  the  signal-to-noise  ratio  (SNR)  that  results  from  a  given  set  of  radar 
and  pulse  parameters.  Knowledge  of  the  SNR  is  necessary  to  allow  us  to  add  realistic 
noise  when  modelling  a  returned  signal. 

Consider  a  coherent  processing  interval  made  of  N  emitted  pulses.  These  pulses  are 
effectively  combined,  or  integrated  coherently,  by  the  operation  of  the  Fourier  transform 
of  Sections  5  and  6  that  produces  Doppler  information  from  those  pulses. 

The  central  equation  that  governs  the  combination  of  the  pulses  is  the  DFT  (6.1). 
When  N  is,  say,  32,  the  DFT  is  a  weighted  sum  of  the  samples  x0, . . .  ,x%i,  and  returns 
the  numbers  X_16, . . . ,  X15.  The  weights  are  unit-magnitude  complex  numbers.  When  N 
complex  noise  samples  are  transformed  in  this  way,  their  weighted  sum  can  be  treated  as 
a  “random  walk”.  It’s  well  known  that  the  root-mean-square  distance  from  start  to  end  of 
N  unit-length  steps  in  a  random  walk  equals  y/N'.  Now  use  the  fact  that  the  power  of  an 
electromagnetic  wave  is  proportional  to  the  square  of  its  amplitude,  so: 

mean  power  of  N  noise  samples  oc  mean  of  (displacement  due  to  N  steps)2 

=  (root-mean-squared  displacement  due  to  N  steps)2 

=  (Vn'  x  distance  covered  by  1  step)2 

oc  N  x  power  of  1  noise  sample.  (El) 

The  two  proportionality  constants  above  are  reciprocals  of  each  other,  so  we  arrive  at 

mean  power  of  N  noise  samples  =  N  x  power  of  1  noise  sample.  (E2) 

This  might  be  expected:  it  says  that  10  light  bulbs  combine  to  give  10  times  the  brightness 
of  one  light  bulb,  because  the  non-coherent  light  from  an  incandescent  bulb  can  be  treated 
as  noise. 

Adding  signals  instead  of  noise  is  a  different  matter,  because  the  phasors  that  describe 
the  signals  are  coherent:  adding  them  is  essentially  like  constructing  a  path  from  the  steps 
of  a  sober  man,  not  a  drunk  one.  The  phasors,  like  the  sober  man’s  steps,  all  point  in  about 
the  same  direction,  so  that  the  distance  from  start  to  end  of  N  of  these  signal  phasors 
is  N  times  the  length  of  one  signal  phasor.  That  means  the  power  present  in  the  sum  of 
N  signals  is  N 2  times  the  power  of  one  signal.  In  other  words,  10  “phase-locked”  lasers 
add  their  light  to  produce  a  central  spot  whose  brightness  is  100  times  the  brightness  of 
the  spot  due  to  a  single  laser  -albeit  that  bright  spot  is  narrower  than  the  spot  due  to  a 
single  laser,  as  expected  from  energy  conservation. 

The  purpose  of  the  weights  e_i 27rkn/N  in  the  DFT  (6.1)  is  to  cancel  the  phase  increases 
that  result  from  any  particular  frequency.  For  example,  if  we  wish  to  extract  from  a 
signal  the  sinusoidal  component  whose  phase  is  increasing  (its  phasor  is  rotating)  at  one 
radian  per  second,  then  we  can  subtract  from  each  sample  a  phase  that  grows  by  one 
radian  per  second,  and  add  the  results.  Only  the  sought-after  component  of  the  signal  will 
then  yield  a  series  of  complex  numbers  of  constant  phase,  so  that  they  add  coherently;  all 
other  components  will  yield  a  series  of  complex  numbers  of  varying  phases,  and  these  will 
produce  some  cancellation  when  they  add.  The  result  is  that  the  sought-after  component 
is  amplified,  while  other  components  are  suppressed.  In  this  way,  the  DFT  adds  signals 
and  noise.  The  peaks  that  result  in  the  range-Doppler  plot  have  their  signal  increased  by 
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a  factor  of  N 2  and  their  noise  increased  by  a  factor  of  N.  So  the  signal-to- noise  ratio  of 
N  pulses  is  N2  /N,  or  N,  times  the  SNR  of  one  pulse. 

The  term  coherent  pulse  integration  denotes  this  IV- fold  gain  from  the  Fourier  trans¬ 
form  that  generates  Doppler  information.  The  signal-to-noise  ratio  of  the  coherently  inte¬ 
grated  power  in  N  pulses  is 


SNR  of  N  pulses  =  N  x 


signal  power  in  one  pulse 
noise  power  in  one  pulse 


We  need  the  signal  and  noise  powers  in  one  pulse. 


(E3) 


Signal  power  in  one  pulse: 

energy  received  in  one  pulse 

signal  power  m  one  pulse  = - — — - - - - —  .  (E4) 

width  of  correlation  peak 

The  radar  has  gain  G ,  carrier  wavelength  A,  losses  L,  and  the  scatterer  has  cross  section  a 
and  distance  r.  The  two-way  radar  equation  says 


energy  received  in  one  pulse 


energy  emitted 
in  one  pulse 


G2aX2 

(47r)3r4L 


(E5) 


The  width  of  the  correlation  peak  is  1/B  where  B  is  the  emitted  pulse’s  bandwidth. 


Noise  power  in  one  pulse:  This  is  kTsB,  where  k  is  Boltzmann’s  constant,  Ts  is  the 
radar’s  system  temperature  (which  is  its  noise  factor  times  the  actual  temperature),  and 
B  is  the  receiver  bandwidth  [6] ,  which  equals  the  emitted  pulse’s  bandwidth  B  if  there  is 
no  reason  to  make  it  smaller  or  larger  than  B:  smaller  would  prevent  the  receiver  from 
receiving  all  that  it’s  required  to,  and  larger  would  make  the  receiver  produce  excess  noise. 
Equation  (E3)  becomes 


SNR  of  N  pulses  =  N  x 


=  N  x 


energy  received  in  one  pulse 
width  of  correlation  peak  x  kTsB 


energy  emitted 
in  one  pulse 


G2aX2 


(4ir)3r4L±kTsB 


energy  emitted 
in  one  CPI 


G2a  X2 

(47t  )3r4L  kTs 


(E6) 


The  gain  of  a  basic  “one-lobe”  radar  is 


G 


47 r 

solid  angle  of  beam 


(E7) 


The  solid  angle  of,  say,  a  beam  of  rectangular  cross  section  just  equals  the  product  of  its 
two  defining  angles,  each  in  radians. 
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E.l  Simulating  Receiver  Noise 


Equation  (E6)  can  be  used  to  derive  an  expression  for  the  noise  which,  when  added  to  a 
pulse,  simulates  the  noise  generated  in  a  receiver. 

First,  consider  two  sources  injected  into  a  narrow-band  IF  amplifier:  a  zero-mean 
gaussian  noise  voltage  of  standard  deviation  u  (we  have  already  used  “<r”  for  cross  section), 
and  a  sinusoidal  signal  voltage  of  amplitude  A.  The  signal-to-noise  ratio  for  one  pulse  is 

mean  signal  power  mean  of  square  of  signal  voltage  A2 /2 

mean  noise  power  mean  of  square  of  noise  voltage  v2 

Referring  to  (E6)  for  N  =  1,  we  then  have,  for  a  transmitted  power  Pt  (considered  constant 
over  one  pulse  width  r,  not  one  PRI), 


A_ 

2v 2 


SNR  for  one  pulse  =  Ptr  x 


G2a\2 

(47t)3  r4 L  kTs  ’ 


(E9) 


from  which  it  follows  that 

Ar2  I  (47 r)3LkTs  ' 
V  =  GAV  2  Ptro 


(E10) 


From  this  v  a  complex-number  noise  can  now  be  generated:  it  has  an  amplitude  drawn  from 
the  gaussian  distribution  AT(0,v2),  and  a  random  phase  uniformly  distributed  between  0 
and  27 r. 


E.2  Maximum  Detectable  Range  of  a  Target 


Suppose  we  know  the  signal-to-noise  ratio  SNRgiven  that  is  expected  for  our  radar,  based 
on  standard  theory  in  radar  texts  that  combines  the  probability  of  false  alarm,  probability 
of  detection,  and  model  of  the  background  clutter  such  as  a  Swerling  model.  Equation  (E6) 
can  use  SNRgiven  to  calculate  the  maximum  range  at  which  our  radar  will  detect  a  target. 
To  see  how,  note  that  in  (E6)  as  the  target-radar  separation  r  reduces,  the  SNR  increases, 
which  is  quite  reasonable.  If  we  now  increase  r  again,  the  SNR  will  drop  until  it  reaches 
SNRgiven.  This  value  of  r  is  the  maximum  range  rmax  at  which  the  target  can  be  detected 
for  the  given  signal-to-noise  ratio  SNRgiven.  We  require  to  find  r,  given  SNRgiven.  Call  the 
energy  emitted  in  one  CPI  “ECPI”,  so  that  (E6)  becomes 


SNR, 


given 


ECP1  G2a\2 
(47r)3r^axL  kTs  ' 


(Ell) 


The  losses  L  might  be  a  function  of  range,  so  rearrange  (Ell)  to  solve  for  rmax: 


Mrn 


ECpi  G~cr \- 


(47t)3  SNRgiven  kTs 


(E12) 


As  an  example,  we  calculate  the  maximum  range  at  which  a  target  of  cross  section 
a  =  5000  m2  can  be  seen  by  a  radar  with  the  following  parameters: 


SNRgiven  =  63  ,  carrier  =  9.3  GHz  (i.e.  A  =  0.0322  m) ,  Pt  =  300  W , 
r  =  16/rs,  16  pulses/CPI ,  G  =  500 ,  Ts  =  710  K ,  PRI  =  500  /rs.  (E13) 
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Note  that  a  power  of  Pt  =  300  W  means  the  300  W  is  treated  as  constant  over  r,  the  pulse 
width;  Pt  is  not  the  mean  power  over  a  PRI. 

For  the  losses  L ,  we’ll  include  only  atmospheric  attenuation.  A  textbook  figure  for  this 
is  0.022  dB/km  at  10  GHz,  which  is  close  enough  to  the  carrier  for  our  calculation.  We 
must  be  very  careful  with  the  units.  I  think  it’s  a  very  under-utilised  observation  about 
units  that  just  as  “6/2”  means  “the  number  of  2s  in  6”,  or  3,  the  expression  “2r/(l  km)” 
means  “the  number  of  kilometres  in  2 r”,  which  is  another  way  of  saying  “2r  expressed  in 
kilometres”.  This  observation  allows  a  sentence  about  units  to  be  expressed  in  concrete 
mathematics.  First,  convert  the  decibels  here  to  bels  to  write  the  attenuation  as  0.0022 
bels  per  kilometre.  This  figure  then  allows  (actually  defines )  the  attenuation  over  our 
two-way  radar  trip  of  distance  of  2 r  to  be  written  as 

L(r)  =  10°'0022T^m.  (E14) 


Equation  (E12)  becomes 


this  is  E, 


CPI 


0.0022 2”max  300  J/s  x  16/rs/pulse  x  16  pulses/CPI  x  5002  x  5000  m2  x  0.03222  m2 


r  10 

'  max 


(4vr)3  x  63  x  1.38  x  10"23  J/K  x  710  K 


~  8.13  x  1019  m4. 

Take  the  fourth  root  of  each  side  to  give 


(E15) 


100°011ft^  ~  95  km, 


(E16) 


and  then  divide  both  sides  by  1  km: 

^lO°-0011^~95.  (E17) 

1  km 

As  expected,  this  equation  is  dimensionless:  we  can  set  a  dimensionless  x  =  rmax/(l  km) 
(“x  is  rmax  in  kilometres”)  and  solve  x  lo°-0011:E  =  95.  This  equation  is  easy  to  solve  by  trial 
and  error;  however,  for  a  sophisticated  approach,  try  using  the  Method  of  False  Position. 
More  generally,  False  Position  never  fails,  unlike  the  more  well-known  but  potentially 
chaotic  Newton-Raphson  method,  which  can  fail  catastrophically  in  general  use. 

The  result  is  x  —  78:  the  maximum  detectable  range  of  the  target  is  rmax  —  78  km. 
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Appendix  F  How  to  Amalgamate  the  Cross 
Sections  of  Many  Scatterers 

The  previous  appendix’s  calculation  of  the  maximum  range  at  which  a  ship  can  be  seen 
required  knowledge  of  the  ship’s  cross  section  a.  In  practice  we  might  know  the  cross 
sections  of  a  thousand  individual  scatterers  that  are  used  to  model  the  ship.  How  can  we 
amalgamate  these  to  form  a  good  estimate  of  the  total  cross  section  a  of  the  ship? 

Refer  to  page  12,  and  consider  the  electric  field  of  the  returned  wave  that  is  incident 
on  the  radar.  As  usual,  we  assume  that  a  single  polarisation  is  being  used.  If  the  electric 
field’s  amplitude  and  phase  are  represented  by  the  complex  number  A,  then  (2.29)  and 
the  comment  after  it  say  that  the  average  areal  power  density  is  (S)  =  e0c|A|2/2.  This  is 
the  average  power  per  unit  area  incident  on  the  receiver.  The  radar  equation  tells  us  that 
this  power  density  is  proportional  to  the  target  cross  section,  so  write 

|A|2  =  K2a  (FI) 


for  some  positive  constant  K  that  depends  on  the  parameters  present  in  the  radar  equation 
(which  is  easily  calculated  but  which  we  don’t  need  to  know). 

But  the  same  argument  applies  to  the  power  density  received  from  the  individual 
scatterers  that  make  up  the  target.  Writing  the  electric  field  amplitude  and  phase  from 
the  jth  scatterer  as  complex  Aj,  we  have 

A  =  J2Aj,  and  \Aj\2  =  r2(Jj  ■  (F2) 

3 


This  last  equation  allows  us  to  write  Aj  =  K for  some  phase  4>j.  Combining  (FI) 
and  (F2)  gives 


kv  =  w2  =  | 

3  3 

The  total  effective  cross  section  of  the  target  is  then 


(F3) 


cr  = 


E 


/a3  e 


Ml 


3  k 


=  E  Va3ak  ei(<^  M 

jk 


E  a3  +  E  \Ja3ak 

3  jk 

j¥=k 


gi {<t>j—<t>k) ' 


(F4) 


For  a  target  composed  of  many  scatterers,  the  phases  —  q vary  randomly  over  all 
values  0  to  27t.  It  follows  that  if  the  scatterers  all  have  similar  cross  sections,  the  last  sum 
over  j  and  k  is  approximately  zero.  Hence  to  a  high  approximation, 


<7  ~  E  a3  ■  (F5) 

3 

That  is,  the  total  cross  section  of  the  target  is  effectively  just  the  sum  of  the  cross  sections 
of  the  individual  scatterers  that  make  it  up.  This  is  a  simple  result,  but  it  requires  a  large 
number  of  scatterers  of  roughly  similar  cross  sections  that  are  all  returning  signals  with 
uncorrelated  phases. 
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Appendix  G  Examples  of  Calculating  Ambiguity 

Functions 

The  ambiguity  function  introduced  in  Section  4.1  is  nothing  more  nor  less  than  the  corre¬ 
lation  of  an  emitted  signal  with  its  return,  when  that  return  has  acquired  some  Doppler 
shift.  The  function  can  be  considered  as  a  stacking  of  correlation  plots,  one  for  each  value 
of  the  target  Doppler,  where  this  stacking  creates  a  surface  in  three  dimensions. 

Specifically,  each  slice  of  this  surface  at  some  target  Doppler  uD  is  the  absolute  value 
of  the  correlation  of  the  signal  with  a  copy  of  itself  that  has  been  Doppler  shifted — 
meaning  each  of  the  copy’s  elements  has  been  multiplied  by  a  Doppler  factor  e^Dt,  where 
t  increments  from  one  element  to  the  next  within  the  copy.  For  example,  the  zero-Doppler 
slice  of  the  ambiguity  function  is  the  absolute  value  of  the  correlation  of  the  signal  with 
an  unaltered  copy  of  itself;  i.e. ,  the  slice  is  the  signal’s  autocorrelation. 

In  principle,  the  signals  that  will  most  clearly  show  the  target  are  those  whose  ambi¬ 
guity  functions  have  a  high  central  peak  and  low  side  lobes.  For  example  in  the  next  few 
pages,  examining  whatever  central  peak  may  be  present,  and  whether  it’s  surrounded  by 
other  peaks,  shows  why  rectangular  non-chirped  pulses  are  not  as  useful  for  radar  signals 
as  are  more  structured  pulses,  such  as  those  that  have  been  Barker  coded. 

Each  of  the  following  plots  shows  the  ambiguity  function  of  the  chosen  signal  type.  For 
each  frequencies  in  a  set  of  Doppler  frequencies,  we  plot  the  absolute  value  of  correlation 
versus  sample  number  using  the  ideas  of  Section  4.2,  and  then  we  stack  each  of  these  plots 
together  along  the  frequency  axis.  This  produces  a  surface  plot  which  is  the  ambiguity 
function  for  the  chosen  waveform.  The  correlation  is  that  of  the  pulse  with  a  copy  of  itself 
that  has  been  Doppler  shifted. 

The  kernel  of  creating  these  ambiguity  functions  is  the  following  Matlab  code,  shown 
here  for  the  case  of  a  single  rectangular  pulse  (although  the  actual  code  used  for  the 
following  plots  contained  more  detail): 

t  =  0  :  0.1  :  4; 

signal  =  ones (1 . length(t)) : 

dopplerFrequencies  =  -5  :  0.1  :  5; 

7,  Initialise. 

correlations  =  zeros (length (dopplerFrequencies) ,  2*length.(signal)-l) ; 

rowNumber  =  0; 

f  or  angularDoppler  =  2*jo_i  *  dopplerFrequencies 
rowNumber  =  rowNumber  +  1; 
cor re  1  at i on ( rowNumber ,: )  =  ... 

Correlate(signal ,  exp (li*angularDoppler*t) ,*signal) ; 

end 

The  matrix  abs (correlation)  can  now  be  plotted.  The  above  code  calls  on  the  one- 
line  function  Correlate,  which  gives  its  output  in  the  order  used  in  this  report,  being  the 
conjugated  reverse  of  Matlab’s  xcorr  output  without  the  unnecessary  extra  zeroes  present 
in  the  output  of  xcorr: 
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function  correlation  =  Correlate (arrayl ,  array2) 

7,  Correlates  arrayl  with  array2  . 

correlation  =  conv(fliplr(conj (arrayl)) ,  array2); 

Figure  G1  plots  the  ambiguity  functions  of  various  pulse  types.  In  particular,  the 
frequency  of  the  chirp  has  been  made  to  change  very  quickly  throughout  the  pulse  to 
highlight  the  rotation  of  the  main  lobe  of  the  surface,  as  compared  to  the  lobe  of  the 
single  rectangular  pulse.  The  Barker  pulse  has  a  well-defined  central  peak,  but  at  a  cost  of 
added  range  side  lobes  that  will,  of  course,  add  some  noisy  structure  to  a  range-Doppler 
plot. 


Rectangular  pulse 


5  rectangular  pulses 


Barker-4 


Chirp 


Figure  Gl:  Ambiguity  functions  for  various  pulse  types:  single  rectangular 
pulse,  a  series  of  five  rectangular  pulses  (whose  PRI  is  twice  the  pulse  width), 
a  Barker-4  pulse  and  a  fast-changing  chirp. 
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repetition  rate,  6 

sampling  interval  restrictions,  42 
speed  measurement  bounds,  33 
speed  resolution,  33 
super-het  receivers,  4 

windowing  data,  27 

z-transform,  67 

zero  padding  in  FFT,  24,  73 


UNCLASSIFIED 


85 


DSTO-TN-1386 


UNCLASSIFIED 


This  page  is  intentionally  blank. 


86 


UNCLASSIFIED 


Page  classification:  UNCLASSIFIED 


DEFENCE  SCIENCE  AND  TECHNOLOGY  ORGANISATION 
DOCUMENT  CONTROL  DATA 


1.  CAVE AT/PRI VAC Y  MARKING 


2.  TITLE 

3.  SECURITY  CLASSIFICATION 

How  to  Create  and  Manipulate  Radar  Range- 

Document 

(U) 

Doppler  Plots 

Title 

(U) 

Abstract 

(U) 

4.  AUTHOR 

5.  CORPORATE  AUTHOR 

Don  Koks 


Defence  Science  and  Technology  Organisation 
PO  Box  1500 

Edinburgh,  SA  5111,  Australia 


6a.  DSTO  NUMBER 

DSTO-TN-1386 


6b.  AR  NUMBER 

AR-016-185 


6c.  TYPE  OF  REPORT 
Technical  Note 


7.  DOCUMENT  DATE 

December  2014 


8.  FILE  NUMBER 

9.  TASK  NUMBER 

10.  TASK  SPONSOR 

11.  No.  OF  PAGES 

12.  No.  OF  REFS 

2012/1057198/1 

07/139 

DMO/MEWSPO 

85 

6 

13.  URL  OF  ELECTRONIC  VERSION 

http : //www . dsto . defence . gov . au/ 

publications/scientif ic .php 


14.  RELEASE  AUTHORITY 

Chief,  Cyber  &  Electronic  Warfare  Division 


15.  SECONDARY  RELEASE  STATEMENT  OF  THIS  DOCUMENT 
Approved,  for  Public  Release 


OVERSEAS  ENQUIRIES  OUTSIDE  STATED  LIMITATIONS  SHOULD  BE  REFERRED  THROUGH  DOCUMENT  EXCHANGE,  PO  BOX  1500, 
EDINBURGH,  SA  5111 


16.  DELIBERATE  ANNOUNCEMENT 

No  Limitations 


17.  CITATION  IN  OTHER  DOCUMENTS 

No  Limitations 


18.  DSTO  RESEARCH  LIBRARY  THESAURUS 


Radar  jamming 
Electronic  warfare 
Signal  processing 


Electronic  countermeasures 
Cross  correlation 
Fourier  analysis 


19.  ABSTRACT 

This  report  lays  out  the  mathematical  framework  and  reasoning  involved  in  addressing  the  question 
of  how  to  produce  false  targets  in  both  range  and  Doppler  to  jam  a  radar.  The  explanations  here 
are  given  from  first  principles  in  order  to  clarify  the  approach  used  in  writing  software  to  address  the 
task.  They  constitute  standard  radar  theory  as  seen  through  the  eyes  of  a  physicist,  and  are  recorded  as 
background  to  the  approach  taken  in  addressing  the  jamming  problem  in  a  later  report,  to  be  published. 

We  discuss  how  a  radar  generates  a  range-Doppler  plot,  using  a  set  of  parameters  that  describe  the 
outgoing  radar  signal  and  the  relevant  characteristics  of  the  target  whose  impulse  response  is  known.  We 
also  explain  necessary  concepts  from  a  mathematical  physicist’s  viewpoint:  bounds  on  pulse  parameters, 
correlation/convolution,  the  theory  of  Hilbert  transforms,  the  relevant  Fourier  analysis,  and  general 
concepts  of  radar  signal  processing  such  as  ambiguity  functions  and  the  maximum  detectable  range  of 
a  target.  This  entire  discussion  is  aimed  at  indicating  the  approach  and  philosophy  used  to  solve  the 
various  problems  encountered  while  working  on  the  task. 


Page  classification:  UNCLASSIFIED 


